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Introduction 

The broad, long-term objective of this IDEA proposal is to achieve optimized image 
quality for DM within acceptable limits of radiation exposure by developing innovative 
approaches for controlling DM exposures. These approaches entail using the digital 
detector and an artificial neural network to control mammographic exposures. This 
project's specific aims are (1) to use short, low dose pre-exposures of the breast to create 
“intelligent” regions of interest that determine the exposure parameters for the fully 
exposed image; and (2) to use an artificial neural network to select exposure parameters 
(mAs, kVp, and beam filtration) based on “intelligent decisions” that optimize signal-to- 
noise (SNR) as a function of mean glandular dose. 


Body: Progress on Research Tasks & Key Research Accomplishments 

Task 1. Determine optimal pixel binning factor and size of exposure-controlling 
ROI. 

The ideal exposure-controlling ROI (ECROI) is large enough so that small, high-density 
structures such as calcifications are not selected to determine exposure parameters, but 
small enough so that the image is not smoothed excessively. To determine an optimal 
ECROI size, we surveyed FFDM images from our clinical image database, evaluating 3- 
D surface plots to determine where the “peaks” of low exposure areas are located and the 
range of the super-pixel values in the peaks. This provided information on the effects of 
the location and size of the ECROI. We evaluated the effect of the ECROI size on the 
calculation of the SNR (Figures 1-7). Next, we developed software to analyze selected 
ROIs (super-pixels) of digitally acquired mammograms. The primary two questions to be 
addressed by these analysis tools are: 

1) Which area of the digitally obtained mammogram contains the area of greatest 
radiographic attenuation (the ECROI)? 

2) What portion of the pixel variance in the ECROI is a result of differential 
breast attenuation (signal) as opposed to random fluctuations (noise)? 

With respect to question 1), the immediate objectives of the software programs are: 

a) Evaluate the correlation between the number of regions of interest (ROIs) 
sampled and the correct identification of the ROI containing the area of greatest 
radiographic attenuation (the ECROI) 

b) Evaluate the correlation between ROI size and correct identification of the 
ECROI 

c) Examine the effect of pixel binning on correct identification of the ECROI 

The selection of the pixel binning factor involves tradeoffs between readout speed and 
information loss. We evaluated binning factors of 8 (to form 320p x 320p super-pixels) 
and 16 (giving 640p x 640p super-pixels) using images of the test phantom acquired at 
low exposures. Three-dimensional pixel-value surface plots were compared to similarly 
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binned plots made using full exposures (Figures 1-4). The location of areas having low 
ADU values was the same in the low- and high-exposure plots, indicating that these 
binning factors did not result in unacceptable image smoothing. We then examined ROIs 
of 800 x 800 pixels (3200p X 3200p pixels) (Figures 5-7) and 1600 x 1600 pixels 
(64,000p X 64,000p. super-pixels). For the 3200p X 3200p pixel ROI, a 6 x 8 grid of 
ROIs was established. In our experiments (Figures 5 & 6), the ROI positions were held 
fixed while their side lengths varied between minimum values of approximately 1.7% of 
the short dimension of the image, to their maximum value, equal to the center-to-center 
ROI spacing. Average pixel value in the 48 ROIs (Figure 5) were plotted as a function of 
ROI side length, in units of binned pixels (Figure 6). We then performed similar 
experiments for 64K p X 64K p ROIs. These data were then utilized for task 2 
experiments (below) to calculate the contribution of signal to the total variance of the 
ECROIs. 

Task 2. Identify the contribution of signal to the total variance in the ECROI. 

This entails distinguishing between contributions to the variance within the ECROI from 
breast structure and from noise (e.g., x-ray quantum noise or a quantum , and additive/system 
noise or a e ) in the pre-exposure images. We have used two pre-exposures to permit pixel 
fluctuations arising from the signal relating to breast structure to be distinguished from 
those arising from noise. The ECROI is determined from the first pre-exposure. Only 
detector modules containing the ECROI are read out on the second pre-exposure. Pre¬ 
exposure times and initial kVp are scaled according to compressed breast thickness (>6.5 
cm:30 kVp; 6 cm:29 kVp; 5.5 cm:28 kVp; 5 cm:27 kVp; 4.5 cm:26 kVp; <4.0 cm:25 
kVp). Initial choices for pre-exposure times and kVp are based on values calculated to 
result in a MGD of approximately 150 mrads to a 50% fat / 50% fibroglandular breast in 
an exposure time of approximately 1 sec. 

Criteria for optimization of tube voltage and external filtration in digital mammography 
differ from those used in screen-film mammography. This is because the separation of the 
processes of acquisition and display in the former permits contrast of individual 
structures to be adjusted when the image is displayed. It is therefore possible to detect 
objects with low subject contrast provided that the image signal to noise ratio (SNR) is 
adequate. Thus, rather than maximization of contrast within the constraints of acceptable 
film darkening and patient dose, beam optimization in digital mammography requires 
maximization of the image SNR, constrained by acceptable patient dose. 

To identify optimum technique factors, the following figure of merit was used: 

FOM = (SNR) 2 /MGD, 

Where MGD is mean dose to the glandular portion of the breast, and the SNR is Rose’s 
suggested minimum value of 5 [1]. 

Using simple phantoms with known signals (disks, holes, etc.) that result in signal- 
induced variances that are easily determined, we calculated the ECROI SNR. Breast 
tissue equivalent material corresponding to 3 different fat/glandular ratios was used to 
simulate a range of breast thicknesses and densities. Look-up tables were generated with 
starting pre-exposure time and kVp for varying breast composition/thickness (288 test 
conditions) as input data. These data were further refined in task 3 experiments using 
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anthropomorphic breast phantoms that more closely approximate the spectral content of 
human breast tissue (see below). 


Task 3. Determine how breast structure relates to measured variances. We selected 
a broad range of breast thickness and radiodensity from our Digital mammography 
patient image database to empirically determine the contribution to the total variance that 
arises from breast structure (signal). In order to obtain the best possible estimate of the 
variance associated with normal breast parenchyma, we selected ROIs from sub-areas of 
mammograms with obviously high SNR. In those ROIs, the average pixel value and 
system gain was used to calculate the x-ray quantum variance, cjq uantum . Using the 

measured system noise, a e , the contribution to the total variance, a total, from breast 
tissue was calculated from:. 


CT signal - CT total a quantum CT e 

After the 2 pre-pulses, we calculated a no ise and a signal and used the empirically 
determined relationship between cr 2 S jgnai and Imax - Imin to infer the magnitude of the 
signal. The ratio of the inferred signal to a noise is the SNR following the pre-exposure. 
This was then used to select exposure parameters for the full exposure. If only the 
exposure time is to be selected (phototiming mode), this is straightforward since I ma x - 
I m j n is directly proportional to mAs and since the quantum noise, a q , increases as the 

square root of the exposure time. In the CCD-based TREX system, the system noise (a e ) 
increases as the square root of the exposure time because of the accumulation of thermal 
electrons in the CCD pixels. This noise was previously characterized over a large range 
of exposure times. 


Results of Phantom Experiments Performed for Tasks 2 & 3: 

Following are example plots of measured signal, SNR, and the FOM, as a function of 
kVp, for a 70/30 fibroglandular/fat equivalent phantom composition, using the signal 
from the microcalcification-equivalent material. These data were obtained from the DM 
system at Johns Hopkins. The results of these studies were presented at the 5 th 
International Workshop on Digital Mammography, Toronto, CA, June 11-14, 2000. 
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Nine different phantoms were assembled to simulate breasts of three different thicknesses 
(3 cm, 5 cm, and 7 cm), and three different attenuation equivalent adipose/fibroglandular 
mass ratios (30/70, 50/50, and 70/30). All blocks of a given phantom had the same 
adipose /fibroglandular ratio, except for two 5 mm thick blocks, common to all phantoms, 
that are 100% adipose equivalent. These blocks were placed at the top and bottom of the 
stack to simulate skin (see figure 4). In each phantom stack assembled, the centrally 
located block in the stack (the signal block) contained a series of test objects. For the data 
reported here, the test objects of interest were two stepwedges, one each of calcification 
equivalent and mass equivalent material. The mass equivalent stepwedge has the same x- 
ray attenuation as 100% glandular equivalent material, and the microcalcification 
equivalent step wedge is composed of calcium carbonate. Figure 5 is a schematic of a 
signal block showing the dimensions of the block and step wedges (other test objects 
present in the signal block have been omitted for clarity). The thickness of all signal 
blocks is 2 cm. Images were obtained in manual mode with the phantoms positioned at 


7 
























































FINAL Report: DAMD 17-99-1-9429 
L.L.Fajardo, M.D. 


the chest wall edge of the receptor, centered left to right. Exposure time was selected to 
give approximately the same average pixel value in the phantom background area for 
each phantom/technique combination. For each combination two images were obtained 
with identical exposure times for the purpose of image subtraction, taking care not to 
move the phantom between the two exposures. At each site, entrance exposures 
(mR/mAs) and half value layers (HVLs) were measured for each target/filter/kVp 
combination used. 



Figure 4: Side view of a 5 cm thick phantom, comprised of one 2 cm thick signal block, 
two 1 cm thick blank blocks, and two 0.5 cm thick skins. 
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Figure 5: Schematic diagram of a signal block 


ECROI Analysis: 

Signal was defined as the difference between the average pixel values in a region of 
interest (ROI) centered on an individual step (but not including the step boundaries), and 
an equal sized ROI located immediately adjacent to the step, but containing only 
background. To quantify the image noise, the two images of a given phantom, obtained at 
a common technique, were subtracted. Image subtraction was performed to remove fixed 
pattern noise associated with phantom defects, detector non-uniformity, and the heel 
effect. Noise in a single image was defined as the rms pixel-to-pixel fluctuations in an 
ROI of 1109 x 511 pixels in the difference image, divided by the square root of two. 
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Evaluation of acceptable mean glandular dose: 

The MGD for each phantom was calculated using its known thickness and composition, 
and the measured HVL and mR/mAs values from each DM system. For Mo/Mo and 
Mo/Rh spectra, the parameterized dose tables of Sobol and Wu were utilized to obtain the 
glandular dose per unit exposure [2]. For the W/Al spectra, normalized (to entrance 
exposure) MGD values were obtained from the data of Stanton et al. [3]. Their data were 
extrapolated to 3 cm breast thickness, and interpolation between their published HVL 
curves was used to obtain correction factors for the particular glandular volume fractions 
(0.22, 0.40, and 0.61, corresponding to glandular mass fractions of 0.30, 0.50, and 0.70, 
respectively) used in our study. For the W/Rh spectra, the calculations of Boone were 
utilized, interpolating between his published HVL and adipose/fibroglandular 
composition values [4]. All FOM values were obtained by dividing the square of the SNR 
by the MGD, expressed in units of 10' 5 Gy (1 mrad). 

The measured HVL values for the seven specific target/filter combinations tested at the 
three sites, as a function of kVp, are shown in Figure 3. Figure 4 shows the corresponding 
normalized MGD, D g N, calculated for each of the seven spectra, plotted versus the 
measured HVL. Similarly, Figure 8 shows D g N for each target/filter combination tested, 
plotted versus kVp. The general tradeoff between loss of contrast and reduction in MGD 
with increasing kVp is illustrated in Figure 9. In this example, the measured contrast of 
the 0.3 mm thick (thickest) calcification step in shown for the 5 cm thick, 50/50 phantom. 

SNR versus kVp, and the corresponding FOM values vs. kVp have been determined. 
Figures 10 and 11 show the results obtained for the 300 micron thick step of the 
calcification stepwedge in the three 50/50 composition phantoms. To illustrate the 
applicability of these data to objects, the dependence of the FOM on the step thickness 
for both types of stepwedges is presented in Figures 12 and 13. These data are from 
images obtained using a Mo/Mo target/filter combination to image a 5 cm thick, 50/50 
composition phantom. Finally, Figures 14-16 illustrate the effect on the FOM of changing 
breast composition, holding breast thickness fixed. Signals were calculated using the 10 
mm thick mass equivalent step. 



Figure 6: Measured HVLs , plotted versus 
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Figure 7: Normalized mean glandular 
dose versus HVL, assuming a 5 cm thick, 
50/50phantom 


Normalized MGD vs kVp : 5 cm,50/50 



Figure 8: Normalized mean glandular 
dose vs. kVp assuming a 5 cm thick, 
50/50 phantom 


Dose and Contrast vs kVp 



Figure 9: Dose and contrast versus kVp 
using the 0.3 mm calcification step in a 5 
cm thick, 50/50 phantom 
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Figure 10: SNR vs. kVp. 


FOM vs kVp : 50/50 



Figure 11: FOM vs. kVp. 



stepwedge, normalized by the average value for each 
step. The average FOM values ranged from 0.2 (step 
0) to 0.01 l(step 4). Imaging data from the 5 cm 50/50 
phantom. 
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Figure 13: FOM values for the four steps 
of the calcification stepwedge, 
normalized by the average value for each 
step. The average FOM values ranged 
from 1.4 (step 0) to 0.64 (step 3). Data 
are imaging the 5 cm 50/50 phantom. 



Figure 14: FOM vs kVp for 3 cm thick 
phantoms of three compositions. 
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Figure 15: FOM vs kVp for 5 cm thick 
phantoms of three compositions 
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Figure 16: FOM vs kVp for 7 cm thick 
phantoms of three compositions. 


Task 4 . Optimize mAs, kVp, and Filter Selection Based on the Results of the Pre¬ 
exposures. 

Optimization is based on maximization of the figure of merit: FOM=(SNR 2 /MGD). We 
have used the added information provided by the spatially varying signal levels in the 
pre-exposure ROIs to calculate a more refined estimate of the MGD than is possible 
using simply compressed thickness. To do this, the mean pixel value in each ROI is used 
to determine the relative transmission through the breast at that location. In regions of 
uniform thickness (all regions whose entrance surface is in direct contact with the 
compression paddle), variations in x-ray transmission are due to variations in breast 
composition. The region of the breast above each ROI is assumed to consist of two skin 
and subcutaneous layers with a uniform mixture of adipose and fibro-glandular tissue in 
between. The x-ray transmission through each region determines an adipose/ 
fibroglandular composition ratio. Thus, for a given kVp, filter, and mAs, the MGD for 
the region of the breast above each ROI can be calculated separately, and the result 
summed to obtain an average MGD. 

Preliminary Results: Task 4: 

Software has been developed for analysis of digitally acquired mammograms. The 
primary two questions to be addressed by these analysis tools may be generalized as 
follows: 

1) Which area of the digitally obtained mammogram contains the area of greatest 
radiographic attenuation (the ECROI)? 

2) What portion of the pixel variance in the ECROI is a result of differential breast 
attenuation (signal) as opposed to random fluctuations (noise)? 

With respect to question 1), the immediate objectives of these programs are: 

1) Evaluate the correlation between the number of regions of interest (ROIs) 
sampled and the correct identification of the ROI containing the area of greatest 
radiographic attenuation (the ECROI) 

2) Evaluate the correlation between ROI size and correct identification of the 
ECROI 
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3) Examine the effect of pixel binning on correct identification of the ECROI 

Figure 17 below shows an example of a digitally obtained mammogram in which a 6 x 8 
grid of ROIs has been established. The figure shows the mammogram, labeled with 
numbers whose locations correspond to the centers of the ROIs. 



7 8 9 10 11 12 


1 2 5 4 -fog s 

Figure 17: A digitally acquired mammogram, 
annotated with the locations of a 6 x 8 square grid 
of ROIs. The numerals are located at the centers 
of the square ROIs. The ROI positions are held 
fixed, while their side lengths are varied between 
minimum values of approximately 1.7% of the 
short dimension of the image, to their maximum 
value, equal to the center-to-center ROI spacing. 



1, plotted as a function of ROI side length, in units of 
binned pixels. For the example shown, the ECROI (that 
with the lowest average pixel value) is the same for all 
ROI sizes tested. Its center is at the position labeled ‘19’ 
in Figure 17. ROIs whose average pixel value is always 
large and constant (top of the graph) are located outside 
the region of the breast. 


Figure 18 is a plot of the mean values of each ROI, plotted as a function of ROI size. 
The ROI with the lowest value is the ECROI. For the example shown, the ECROI is 
ROI #19. The contour plot of the pixel values, shown below in Figure 19, verifies 
that this is in fact the most highly attenuating portion of the breast. 
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Figure 19: Contour plot of the pixel values in the digital 
mammogram of Figure 1. 
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We have developed analysis tools to evaluate the effect of radiographic heterogeneity on 
the selection of the ECROI and have performed analysis of digital mammograms of 
patients with a variety of compositions and breast sizes. Table 1 below provides the 
calculated ratios of root mean square signal to root mean noise for 20 digital 
mammographic images evaluated in 5 patients with a variety of breast compositions 
(radiodensity). For each study, the compressed breast thickness was 4 -5 cm and a 28 
kVp Mo/Mo beam was used. The third column gives the total variance in the ECROI. 

The fifth column gives the portion of the variance that is due to noise, as determined by 
uniform irradiation studies. The last column gives the ratio of the RMS signal to the RMS 
noise. These ratios suggest the approximate magnitude of the target SNR values to be 
used in the AEC algorithm. 


TABLE 1: Ratio of RMS Signal to RMS Noise 



Average 

Total 

Variance 

mAs 

Noise 

Variance 

Signal Variance 

SNR 

Patient 1 

1422 

17986 

43 

745 

17242 

4.8 


1420 

11556 

43 

744 

10812 

3.8 


1504 

29009 

46 

780 

28229 

6.0 


1630 

23148 

50 

835 

22313 

5.2 

Patient 2 

2031 

36416 

62 

1019 

35397 

5.9 


1964 

16243 

60 

987 

15255 

3.9 


1964 

25290 

60 

987 

24303 

5.0 


1972 

26219 

60 

991 

25228 

5.0 

Patient 3 

1535 

9923 

47 

794 

9129 

3.4 


1430 

16615 

44 

748 

15867 

4.6 


1458 

14053 

44 

760 

13293 

4.2 


1461 

19624 

45 

762 

18863 

5.0 

Patient 4 

1577 

25117 

48 

812 

24305 

5.5 


1502 

20446 

46 

779 

19667 

5.0 


1547 

22630 

47 

799 

21831 

5.2 


1311 

35793 

40 

698 

35095 

7.1 

Patient 5 

2204 

49051 

67 

1102 

47949 

6.6 


1727 

20817 

53 

879 

19938 

4.8 


2041 

22180 

62 

1024 

21156 

4.5 


2329 

39980 

71 

1163 

38816 

5.8 
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We have carefully evaluated linearity as a function of x-ray exposure (Table 2 & Figure 
20a) with phantom experiments using a 4cm acrylic phantom imaged at 28 kVp. The 
quantum limited range of operation was determined from a log-log plot of the pixel 
variance vs the mAs (Figure 20b). As is shown, the linear region between ~ 1.0 and 2.5 on 
the horizontal axis (Figure 20b) indicated quantum limited operation. Below log(mAs) = 

1.0, the system noise constitutes a significant part of the total noise, and above log(mAs) = 
2.5, there is a slight departure from quantum limited operation, probably due to imperfectly 
corrected fixed pattern noise at these high exposure levels. This departure can also be seen 
in the two data points in the upper right portion of Figure 20c. We will use the quantum 
limited region in our calculation of the final mAs, based on the measured signal variance 
(pixel variance due to differential breast attenuation) in the ECROI following pre¬ 
exposure. . 

We have also developed analysis programs for determination of the pixel variance due to 
differential breast attenuation (signal variance), and are characterizing intrinsic detector 
noise. Table 3 below provides data on the total variance in the ECROI, including its 
components due to signal and noise. 


Table 3. Total Variance in the ECROI as a function of the 
Signal and Noise Components 


mAs 

Average 

Total Variance 

Noise Variance 

Signal Var 

sqrt(sig var) 

3 

43 

262 

167 

95 

10 

6 

94 

409 

215 

194 

14 

10 

166 

700 

279 

421 

21 

20 

343 

2110 

448 

1662 

41 

40 

698 

7548 

803 

6745 

82 

80 

1401 

29174 

1507 

27667 

166 

160 

2813 

116440 

2979 

113461 

337 

240 

4224 

261613 

4734 

256879 

507 

325 

5724 

481110 

8617 

472493 

687 

400 

7056 

732557 

11067 

721490 

849 
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FIGURE 21a. Average Pixel Value in ECROI 



Figure 21b. Total Variance in ECROI 



Figure 21c. S.D. in ECROI due to Signal only 

Standard deviation in ECROI due to signal only 
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These results demonstrate that the total variance in the ECROI is the arithmetic sum of the signal 
and noise variances. As shown in Figure 21c, the signal variance is proportional to the square of 
the mAs, as would be expected for a linear system. Therefore, determination of the total variance 
at any given mAs, along with the known noise variance at that mAs, permits calculation of the 
SNR for all mAs values. In the context of AEC development, this means that we can quickly 
assess the RMS signal variance in the ECROI following pre-exposure at low mAs, and 
extrapolate to the target mAs value that will result in the desired SNR. 


Key Research Accomplishments: Progress Since Last Report 

Over the last 12 months, we evaluated radiation exposures and image quality in a population of 
women (n=245) undergoing both conventional film mammography (FM) and DM. Four images 
or “views” were performed by each modality. Information regarding the patients’ age, breast 
tissue radiodensity, the breast thickness when compressed for the mammograms, and exposure 
factors were collected. The radiodensity of breast tissue on a mammogram correlates with the 
ease or difficulty of reading the study; cancers are more easily detected in tissue composed 
entirely of adipose tissue (not radiodense) than in breast tissue that is more radiodense (i.e., 
comprising fibroglandular rather than adipose tissue). 

For each modality, we calculated the mean glandular radiation (MGD) dose for each of the 4 
images comprising each subject’s DM and FM and performed an assessment of image quality 
using radiologists certified in mammography. Here, the traditional FMs using a standard 
radiographic “light box” and DMs were evaluated on a soft copy workstation equipped with high 
resolution CRT monitors for displaying the digital images and with the capability to adjust the 
image brightness, contrast, and window levels for optimal interpretation. For each modality, the 
radiologist readers completed an evaluation form, providing subjective assessments of image 
quality with respect to visualization of critical breast anatomic structures, using a rating scale of 
l=poor to 5=excellent quality. 

Results 

The MGD for each of 4 breast images associated with each case (245 DMs and 245 FMs, 
therefore 980 DM and 980 FM images) was calculated and evaluated by analysis of variance 
with repeated measures for modality and view (Table4). The average MGD per image was 
greater for FM than for DM: 217.1mrad vs. 208.0 mrad (p = 0.0300; F(l,244) = 4.76). The 
magnitude of the MGD difference varied with the particular mammographic view. 


Tal 

ble 4. MGD (mrads) per Mammographic View 

Modality 

Lt MLO 

Rt MLO 

Lt CC 

Rt CC 

FM 

224.8 

227.1 

207.6 

209.1 

DM 

212.8 

212.2 

205.1 

201.9 

p = 0.0214; F(3,732) = 3.54 

MGD=mean glandular dose, mrad=millirads, Lt=left, Rt=right, MLO=mediolateral oblique 
view, CC=craniocaudal view, FM=film mammography, DM=digital mammography 
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The overall MGD (adding together the dose from all four views) for FM and DM was studied by 
analysis of variance with a grouping factor for breast radiodensity using FM as the gold standard 
and a repeated measure for mammographic modality. Overall MGD was higher for FM than for 
DM: 924.5 mrad vs 855.2 mrad (p = 0.0096; F(l,241) = 6.81). The difference in overall dose 
between film and digital DID NOT depend on breast radiodensity (F(3,241) = 1.53, p = 0.2064). 


Overall Mean Glandular Dose 



□ Film 
■ Digital 


^ J 1 

•v ^ «f 

Breast Radiodensity 


Fimire 22 


Table 5. Mean Image Quality 
for FM and DM by Breast Density 

Breast Radiodensity 

FM 

DM 

Not radiodense 

3.9 

4.7 

Minimally radiodense 

3.8 

4.7 

Moderately radiodense 

3.7 

4.7 

Extremely radiodense 

3.0 

4.4 

p = 0.1516; F(3,209) = 1.78 


The image quality ratings (l=poor to 5=excellent) were lower for FM than for DM: 3.6 vs 4.6 (p 

< 0.0001; F( 1,209) = 185.96). The magnitude of the quality advantage for DM depended on the 
specific features being evaluated with the skin line on DM being better visualized than on FM (p 

< 0.0001, F(3,627) = 49.72) (Table 6). However, the advantage for DM did not depend on breast 
radiodensity (Figure 23). 
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Table 6. Mean Quality Rating for FM and DM of Specific Image Features 


Posterior Tissue 

Retroareolar tissue 

Retroareolar Structures 

Skin Line 

Film 

3.6 

3.7 

3.8 

3.2 

Digital 

4.6 

4.6 

4.6 

4.8 
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Reportable Outcomes: 

Abstracts & Presentations: 

The 1999 Radiological Society of North America, Chicago, IL, November 27- December 3, 
1999: “Development of a quality control system for full-field digital mammography”. 

MJ Yaffe, MB Williams, LT Niklason, GE Mawdsley, AD Maidment, Radiology 209(P) (1999) 
160. 

The 5 th International Workshop on Digital Mammography, Toronto, CA, June 11-14, 2000: 
“Beam Optimization for Digital Mammography”. MB Williams, M More, V Venkatakrishnan, L 
Niklason, MJ Yaffe, G Mawdsley, A Bloomquist, A Maidment, D Chakraborty, C Kimme- 
Smith, LL Fajardo. 

Submitted: 2004 RSNA”, November 2004. “Breast Radiation Dose for Conventional Film Versus 
Digital Mammography”, Grignon L, Berbaum K, Williams MB, More MG, Maxwell S, 
McGruder A, Fajardo L. (Included in appendix). 


Publications: 

Jung W, Liu H, Fajardo LL. Contrast detail detectability of a full field digital mammography 
system. SPIE Proceedings 1999;3595:186-192. 

Liu H, Jiang H, Chen W, Fajardo LL. Lens distortion in optically coupled digital x-ray 
imaging. Medical Physics 2000;27(5):906-912. 

Williams MB, Mangiafico PA, Simoni PU. Noise power spectra of images from digital 
mammography detectors. Med Phys 2000:6(7); 1279-1293. 


Conclusions: 

The analysis of SNR and FOM as a function ofkVp, shown in Figures 10 and 11, indicates that 
although the image SNR tends to decrease monotonically for all systems with increasing kVp, 
the accompanying MGD reduction results in fairly flat FOM curves. This is primarily due to tube 
loading, since it was not possible to obtain the same exit exposure at all kVps (that is, the tube 
output was insufficient to compensate for the lower transmission through the phantoms). Thus 
the falling SNR (and the falling MGD) with decreasing kVp are really consequences of falling 
exposure. 

For a given phantom/technique combination, the SNR, and thus the magnitude of the FOM, 
increases with increasing step thickness for both types of stepwedge. However, the shape of the 
FOM vs. kVp curves for a given target/filter/phantom combination are essentially independent of 
step thickness, and are similar for mass and calcification equivalent signals. This is illustrated by 
the example shown in Figures 12 and 13. This implies that the result of the optimization is not 
sensitively dependent on signal amplitude. 
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Figures 14-16 illustrate a clear advantage to using rhodium filtration for thin breasts, but that for 
breasts 5 cm or thicker, aluminum filtration becomes increasingly advantageous. Similar 
statements can be made for the molybdenum target systems tested, where molybdenum filtration 
was superior for 3 cm phantoms of all compositions, but rhodium filtration produced better 
results for 5 and 7 cm thick phantoms of all compositions. These data suggest that the choice of 
external filtration is potentially more significant in determination of the overall FOM of a DM 
system than is choice of tube voltage. 

Fahrig and Yaffe developed a model for optimizing spectral shape in digital mammography, and 
used it to calculate kVp values producing maximum SNR at a fixed dose for W and Mo spectra 
[5]. They found that, for a fixed MGD of 0.6 mGy (60 mrads), the peak SNR occurred in the 24- 
31 kVp range (W spectrum) and 25-29 kVp range (Mo spectrum) for 4 - 8 cm breast thickness, 
and 50/50 breast composition. Their results were the same, whether the lesion type modeled was 
infiltrating ductal carcinoma or microcalcification. 

Jennings et al. used a computational approach to identify maximum FOM values (FOM = 
SNR 2 /MGD) for a variety of target/filter combinations, and breast thicknesses [6]. They found 
that for a Mo/Mo beam used to image 3-6 cm, 50/50 breasts, the FOM peaks at 27-28 kVp, and 
changes slowly with changing kVp near the peak values. Very similar FOM vs. kVp curves were 
obtained for Mo/Mo, Mo/Rh, and W/Al spectra, applied to 6 cm thick, 50/50 composition 
breasts. The general trends in our data appear to be consistent with those of these previous 
studies. 

Our results using the Lorad CCD-based detector demonstrate that the total variance in the 
ECROI is the arithmetic sum of the signal and noise variances, and that knowledge of detector 
noise characteristics for a given attenuating thickness and kVp, permits the RMS variance due to 
differential breast attenuation to be obtained. Because the system is linear with respect to x-ray 
fluence, the signal variance is proportional to the square of the mAs. Therefore, calculation of 
both the signal and noise variance at any given mAs is possible, enabling permitting calculation 
of the SNR as a function of mAs . 

Implementation of the first generation of automatic exposure control device on the digital 
mammography system results in images having satisfactory image quality, SNR and mean 
glandular dose. For 16 patients undergoing both conventional mammography and digital 
mammography, the implemented AEC for the digital system performed satisfactorily. We now 
have nearly complete data to finalize the algorithms for a functional AEC device for the 
experimental digital mammography system. 


Our data indicate a statistically significant trend toward better image quality for DM. Although 
the image quality ratings are subjective, we are encouraged that quality is at least comparable 
and perhaps better than conventional FM. A statistically significant reduction in overall MGD 
and in MGD per view for DM was also measured. Thus, lower radiation exposures can be 
achieved with DM while maintaining subjective image quality. These results encourage further 
investigation into the dose reduction capabilities of DM as this new technology disseminates 
more widely into clinical use for breast cancer screening. 
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Lens distortion in optically coupled digital x-ray imaging 
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The objectives of this research are to analyze geometrical distortions introduced by relay lenses in 
optically coupled digital x-ray imaging systems and to introduce an algorithm to correct such 
distortions. Methods: The radial and tangential errors introduced by a relay lens in digital x-ray 
imaging were experimentally measured, using a lens-coupled CCD (charge coupled device) proto¬ 
type. An algorithm was introduced to correct these distortions. Based on an x-ray image of a 
standard calibration grid, the algorithm first identified the location of the optical axis, then corrected 
the radial and tangential distortions using polynomial transformation technique. Results: Lens dis¬ 
tortions were classified and both radial and tangential distortions introduced by lenses were cor¬ 
rected using polynomial transformation. For the specific lens-CCD prototype investigated, the mean 
positional error caused by the relay lens was reduced by the correction algorithm from about eight 
pixels (0.69 mm) to less than 1.8 pixels (0.15 mm). Our investigation also shows that the fourth 
order of polynomial for the correction algorithm provided the best correction result. Conclusions: 
Lens distortions should be considered in position-dependant, quantitative x-ray imaging and such 
distortions can be minimized in CCD x-ray imaging by appropriate algorithm, as demonstrated in 
this paper. © 2000 American Association of Physicists in Medicine. [S0094-2405(00)03705-6] 

Key words: lens distortion, optically coupled x-ray imaging system, camera calibration 


I. INTRODUCTION 

Both fiber optically coupled and lens-coupled CCD x-ray 
imaging systems have been used in digital mammography 
and digital radiography. 1-3 The cascaded imaging chain con¬ 
sists of a scintillator, an optical component (either a relay 
lens or a spatial coherent fiber bundle), and an electronic 
imager, such as a CCD imager. 4-6 Obviously, the optical 
properties of the lens, or the fiber, can affect the quality of 
the x-ray imaging significantly. In our previous communica¬ 
tions, the characteristics of optical fiber taper and its impact 
on image contrast were analyzed. 7 The impact of lens¬ 
coupling efficiency on signal-to-noise ratio was 
investigated. 8 This communication focuses on the geometri¬ 
cal distortions caused by relay lenses, since it is well known 
that aberrations, including geometrical distortions, do exist 
even with a well-designed and well-built lens. We will first 
introduce the origin of lens geometrical distortion. We will 
then take the lens-coupled CCD x-ray imaging system as an 
example to analyze the positional errors introduced by lens 
geometrical distortion. Finally, we will present an algorithm 
for correcting geometrical distortions caused by lens in digi¬ 
tal x-ray imaging. 

II. GEOMETRICAL DISTORTIONS CAUSED BY 
LENS 

As shown by Fig. 1, a relay lens used in an optically 
coupled CCD x-ray imaging system collects light from a 


point on a scintillator, and focuses it onto a corresponding 
conjugate point on an electronic imager. Under the most 
practical conditions, the lens, same as all other optical com¬ 
ponents, is not able to form a perfect image due to the pres¬ 
ence of aberrations. Lens aberrations include spherical aber¬ 
ration, coma, field curvature, astigmatism, axial color, lateral 
color, and distortion. 9 Most of these aberrations degrade the 
image quality by “blurring” the image. 

Lens distortion is a unique aberration in that it does not 
degrade the quality of the image in terms of sharpness and 
focus. Rather, it affects the shape of the image, causing it to 
depart from a true-scaled duplication of the original object. 
For instance, if one acquires the image of a square grid as 
shown in Fig. 2(a), the recorded image can be affected by 
pincushion distortion, as shown in Fig. 2(b), and affected by 
barrel distortion, as shown in Fig. 2(c). Pincushion is an 
outward displacement of a given image point from its desired 
location on a mean image plane: Its magnitude increases, in 
the simplest cases, as the field angle of the optical system 
increases. The field angle a is defined in Fig. 1. On the other 
hand, barrel distortion is an inward displacement, and its 
magnitude also increases with the field angle, but with a 
different orientation [see Fig. 2(c)]. The distortions demon¬ 
strated in Fig. 2 are only radial distortions, which for a 
purely rotational system cause the actual image point to be 
displaced radially in the image plane, independent of the 
azimuth in the field. 
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Fig. 1. A schematic of a lens-coupled CCD digital x-ray imaging system. 
The field angle a is defined as the angle formed by the optical axis (dotted 
line) and the principal ray through the center of the lens connecting the two 
conjugate points. 


In principle, a lens is a rotationally symmetric system. In 
practice, no ideal lenses can be made. Small errors retained 
in the alignment and centering of a multielement optical sys¬ 
tem lead to tangential displacement of image point; the mag¬ 
nitude of such distortion depends on both the azimuth and 
the off-axis position of the field. Figure 3 shows an example 
of the tangential distortion. One of the effects of tangential 
distortion is that a straight line passing through the center of 
the field may be recorded as a weakly curved line. 

Lens aberrations cannot be avoided completely in optical 
design; however, they can be minimized. Optical designers 
evaluate the potential contribution of each aberration to the 
final system performance and adjust the configuration of the 
optical assembly to achieve satisfactory performance for 
each application. For instance, an optical designer usually 
balances the outward and inward displacements to minimize 
geometrical distortions. The designer, however, has to 
“trade-off’ many factors, such as cost and size, against 
technical parameters required by specific applications. With 
exception of certain custom built systems, distortion errors in 
the range of 1-3 percent may be present in commercially 
available relay lenses, when these lenses are used in the im¬ 
aging system. 



Fig. 2. Illustration of geometrical distortions introduced by lens, (a) A 
checkerboardlike grid without distortion; (b) Pincushion distortion; (c) Bar¬ 
rel distortion. 



Fig. 3. Illustration of tangential distortion: A straight line passing through 
the center (optical axis) of the field is recorded as a weakly curved line. In 
this schematic, the dashed lines show the undistorted geometrical pattern 
and the solid lines represent the distorted pattern. 


Due to the presence of the residual distortion, the local 
scale on an image produced by a lens-coupled digital x-ray 
imaging system may vary both radially and tangentially at a 
given image position from the mean value. 

In many diagnostic-imaging applications, such as routing 
gastrointestinal fluoroscopy, the presence of geometrical dis¬ 
tortions is often unnoticed by the image readers (its effect on 
the diagnostic quality is up to discussion). The presence of 
geometrical distortions is, however, of great importance in 
position-dependent, quantitative imaging, such as in cardiac 
catheter placement, neurological embolization procedures, 
and in image guided biopsy. For instance, the errors intro¬ 
duced by both radial and tangential distortions reduce the 
positional accuracy of a lens-coupled CCD stereotactic sys¬ 
tem. 


III. POSITIONAL ERROR INTRODUCED BY LENS 

To investigate the geometrical distortion introduced by 
lenses in digital x-ray imaging, a lens-coupled CCD proto¬ 
type was used. 3 The system, consisting of a scintillating 
screen (Gd 2 0 2 S:Tb), a relay lens and a CCD camera, is at¬ 
tached to a clinical x-ray machine. The camera employs a 
mechanical shutter. An optical hood was used to shield the 
CCD sensor from ambient light. The CCD array size is 
1024X 1024 pixels, and 0.024 mmX 0.024 mm per pixel. The 
camera was operated under a temperature of —25 °C to re¬ 
duce thermal electron noise. The low temperature was 
achieved with a compact thermoelectric cooler. The opera¬ 
tion of the system is described as follows: (1) During x-ray 
exposure, the x-ray beam passing through the testing target is 
attenuated and interacts with the scintillating screen. (2) The 
x-ray photons are then converted into a large number of vis¬ 
ible light photons. (3) The relay lens projects the optical 
scene from the scintillating screen to the photosensitive sur¬ 
face of the CCD camera. (4) The CCD pixel array samples 
the information and creates a digital image. In our experi- 
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Fig. 4. Digital x-ray images showing the effect of lens distortion, (a) An 
x-ray image of a checker board like grid, acquired using a prototype lens- 
coupled CCD x-ray imaging system, (b) a standard map of the grid (dashed 
lines), properly scaled, was artificially overlaid on the distorted x-ray image, 
to demonstrate the geometrical distortion. 



Fig. 5. Illustration of the coordinate system to characterize lens distortions. 
The optical center (x c ,y c ) is served as the origin of polar domain. The base 
axis of the polar system is parallel to the x axis of Cartesian system. The 
distorted and undistorted points are denoted as (p,0) and (p',0')> respec¬ 
tively. The radial distortion, Ap, tangential distortion, A r, and combined 
distortion. Ay, are also illustrated in this figure. 


ment, this operation was controlled by a computer. The x-ray 
energy and exposure used in the experiments were 30 kVp 
and 32 mAs, respectively. 

A commercial Nikon lens, F/1.2, 50 mm focal length, 
equipped with a close-up (Nikon Close-up #2) lens was used 
as the relay lens to project the optical image from the scin¬ 
tillating screen to the CCD imager. The relay lens is working 
at a 3.6:1 demagnification ratio, thus making an object of one 
millimeter in length occupy — 11.6 pixels in digital image. 
During experiments, a checkerboard-like grid made by cop¬ 
per wires was imaged by the prototype system. To minimize 
the effect of x-ray focal spot on image quality and geometri¬ 
cal shape, the copper wires were kept directly against the 
scintillator (no air gap between the testing target and the 
x-ray detector), and the distance from the x-ray focal spot to 
the scintillator was kept at 650 mm. The diameter of the 
copper wires was 0.17 mm and the spacing between the cop¬ 
per wires (center to center) was 5 mm. 

An x-ray image of the calibration grid acquired by our 
imaging system is shown in Fig. 4(a). In Fig. 4(b), a map of 
the grid (dashed lines), properly scaled, is artificially over¬ 
laid over the acquired image, to demonstrate the geometrical 
distortion. At locations near the center of the image (near the 
optical axis of the lens), the standard map overlaps with the 
x-ray image almost perfectly. At locations toward the edges 
of the image, the displacement between the map and the 
x-ray image becomes clearly noticeable. Apparently, this 
specific lens used in this prototype introduced a significant 
barrel distortion. 

To quantitatively measure positional errors introduced by 
lens distortions in digital x-ray imaging, positional displace¬ 
ment of the grid image from its true location was calculated 
for all the feature points (intersections of the grid wires). In 


Table I, radial distortion errors, and total distortion errors 
(combined radial and tangential errors) are listed as functions 
of the field angle of the optical system. The radial distortion 
error, Ap, and the total distortion error, Ay, are illustrated in 
Fig. 5 in a polar coordinate. 

As shown in Table I, the positional displacement of an 
object imaged by such a lens-coupled digital x-ray imaging 
system can be, in the most severe cases, more than nine 
pixels away from its true location within a field angle of 30°. 
Considering the fact that the CCD imager used in the proto¬ 
type has only 1024 X 1024 pixels, the positional error intro¬ 
duced by the lens distortion cannot be neglected in any 
position-dependent, quantitative x-ray imaging. 


IV. AN ALGORITHM FOR LENS DISTORTION 
CORRECTION 

Digital image processing techniques can be used to cor¬ 
rect geometrical distortions introduced by optical and opto¬ 
electronic components. 10-12 In our study, a correction algo¬ 
rithm was introduced to minimize the lens distortion. The 
algorithm was derived based on the nature of the lens distor¬ 
tion. 

After a digital image of a standard grid is acquired by the 
imaging system, the algorithm operates as follows: 

(1) Localize the intersections of the wires of the grid (fea¬ 
ture points); 

(2) Create an ideal map of the grid; 

(3) Identify the optical center of the lens; 

(4) Perform polynomial-transformation to correct distor¬ 
tions. 


Table I. Positional errors introduced by relay lens versus field angle of the optical system. 


Field angle a (degree) 

3.6 

8.6 

Radial distortion A p (pixel) 

1.26 

2.28 

Combined distortion Ay (pixel) 

2.39 

3.16 

Combined distortion Ay (mm) 

0.20 

0.27 


13.0 

17.4 

22.1 

26.2 

30.0 

3.25 

4.30 

5.43 

6.94 

7.42 

4.15 

5.19 

6.61 

8.04 

9.71 

0.36 

0.45 

0.57 

0.69 

0.84 
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These procedures are discussed in detail in the succeeding 
text. 

A. Localizing feature points 

A semiautomatic procedure was developed to localize the 
positions of the feature points (wire intersections of the grid). 
The first step of this procedure was to manually locate two 
adjacent feature points, which served as the starting points of 
the semiautomatic procedure. 13 Centered at each of the start¬ 
ing points, a local window was set up (40X40 pixels rect¬ 
angle in our experiment). A connected region was deter¬ 
mined that contained m adjacent pixels (m = 8 in our 
experiment) with the lowest gray levels inside the window. 
This region could be named as a valley region. Because of 
the superposition of the two wires, the x-ray attenuation at 
each intersection of the grid (feature point) was higher. 
Therefore, feature points showed minimum pixel value in the 
image. It is reasonable to assume that the feature point was 
located inside the valley region. The reason for chosen a 
connected region instead of a pixel is to improve the algo¬ 
rithm’s robustness for the unavoidable imaging noise, and to 
improve the accuracy of the locations of the feature points. 
The exact location of the feature point ( x a ,y a ) was defined 
as the gravity center of the valley region 

1 m 

*<, = — 2 xJ{Xi,yi) 

m , = i 

1 m 

ya=—^yif(x t ,yi), 0 ) 

m i= i 

where f{x i9 y t ) was the gray level at point (*,,>>,) in the 
valley region. 14,15 Using Eq. (1), the first two feature points, 
their starting points being manually located, can be deter¬ 
mined in sub-pixel accuracy. 

The algorithm could identify the next starting point ac¬ 
cording to the distance between two known feature points 
along the extension (or perpendicular) line between them, 
since the calibration image was a checkerboardlike grid. One 
by one, the remaining feature points were determined by 
using this method. 

B. Creating the ideal map of the grid 

The ideal map of the grid was created to serve as the 
standard pattern for image correction. In Ref. 16, Shah and 
Aggarwal acquired an ideal image by another lens that intro¬ 
duced minimal distortion. If the polynomial warping method 
was employed, the grid target could be used as an undis¬ 
torted ideal grid itself. 16 Here, a method was employed to 
create an ideal map of the grid from the actual distorted grid 
image using the following mechanism. 

The grid pattern consisted of two sets of parallel straight 
lines intersecting at a right angle. One set of parallel lines 
can be described by the equation: 

y t = kx + ib 1,0,1,...), (2) 

where the parameters k and b are the slope and the 
y-separation of the lines, respectively. It is reasonable to as¬ 


sume that the lens distortion is negligible near the center of 
the image. Thus, the parameters k and b can be estimated 
using feature points in the center region of the actual grid 
image. In our experiment, with the whole image size of 
1024X1024 pixels, a 200X 200-pixel sub-image was taken 
around the image center for the purpose of the parameter 
calculation. This sub-image covers four pairs of straight 
lines. Each line equation was initially determined using two 
calculated feature points as described above. Then in the 
perpendicular direction along this line, the points with a local 
minimum gray level were localized. Using the least mean- 
square fitting method, the parameters of the line equation 
could be determined precisely with these local minimum 
gray level points. The calculated grid map of the center re¬ 
gion was then extended to the whole image, thereby creating 
the entire calculated grid map. However, the map was 
roughly positioned with respect to the geometrical center of 
the x-ray image array. The last step of this procedure was, 
therefore, to shift the calculated grid map slightly to achieve 
the best alignment with the actual grid image by minimizing 
the following criterion function 

^=2 [(x g i-x al ) 2 +(y gi -y al ) 2 ], (3) 

i 

where (x g ,y g ) and (x a ,y a ) are the feature points corre¬ 
sponding to the calculated grid map and actual grid image, 
respectively. The subscript i indicates the ith feature point, 
and the summation in Eq. (3) is over all feature points of the 
x-ray image. 

C. Identifying the optical center of the lens 

Normally, the lens distortion is mainly in radial direction, 
which causes an inward or outward displacement of a given 
image point from its true location. This type of distortion is 
strictly symmetric about the optical axis. 17-19 With this prop¬ 
erty, a method was developed to localize the optical center 
(optical axis). 

Let (x c ,y c ) denote the optical center, the radial distortion 
at any feature point (x ai9 y ai ) in the actual image can be 
calculated by: 

d(x a j >yai) V(-^ai X c ) 2 + (y a i~~ y c)^ 

- 'I(x gi -x c ) 2 + (y gt -y c ) 2 , (4) 

where, (x gi9 y gi ) is the corresponding feature point to the 
( x ai ,y ai ) in the calculated ideal grid map. 

For feature point (x ai ,y ai ), its symmetric location about 
the optical center (x c ,y c ) is given by (2x e -x ah 2y c —y ai ). 
Therefore, the asymmetry of the radial distortion to (x c ,y c ) 
can be characterized by a function 

d(x c ,y c ) \.d(x a j 9 y a i) d{2x c x a i,2y c y^.)]* (5) 

i 

where the subscript i indicates the zth feature point. The sum¬ 
mation is over all feature points of the x-ray image. The 
optical center is such a point that the asymmetry of the dis¬ 
tortion is at a minimum level. Therefore, the problem of 
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searching the exact location of the optical center is equiva¬ 
lent to an optimization procedure in which the center posi¬ 
tion (x c ,y c ) is to be determined in order to minimize the 
object function J in Eq. (5). 

To perform such an optimization procedure, the following 
steps were adopted: 

(1) First, select the center of the digital image as the initial 
location for the optical center. 

(2) Search for x c , by minimizing the object function J in 
Eq. (5) with y c fixed. 

(3) Then, with a fixed x c determined in step 2, search for y c 
by minimizing J. 

(4) Repeat steps 2 and 3 until the location of the optical 
center is stabilized within a reasonable tolerance (0.2 
pixels in our experiment). 

It should be noted that for a feature point ( x a ,y a ), its 
symmetric point (2x c —x a ,2y c — y a ) about the center (jc c ,y c ) 
may be falling between the known grid wire intersections. In 
such a situation, a two-dimensional interpolation was imple¬ 
mented according to its neighbors to seek an estimate of its 
radial distortion. Bilinear interpolation was used in our 
method. 20 

D. Performing polynomial-transform to correct 
distortion 

Polynomial transformation is a common method used to 
correct the distortions. 21-23 Since a combination of radial dis¬ 
tortion and tangential distortion was seen in the image, in our 
experiment two sets of polynomials were needed. One was 
used to correct the radius as measured from the optical cen¬ 
ter, and the other was used to correct the polar angle as 
measured from the base axis of the polar system. As shown 
in Fig. 5, the distorted image of a feature point is denoted as 
polar coordinate (p,0), and the distortion-free image of the 
point is denoted as (p',0'). The radial position of each pixel 
and the tangential offset component were corrected by apply¬ 
ing two set of N-order polynomials: 

0'=a l 0+a 2 0 2 + a 3 0 3 + a 4 6> 4 + --- > 

p' = b i p+b 2 p 2 + b i p 2 + b 4 p 4 -\ -, 

where, a t and b t are the distortion coefficients. The above 
equations hold for all feature points of the image. Because 
the number of feature points generally exceeded the number 
of polynomial coefficients, a { and b t were solved by least- 
squares fitting which minimized the mean-square error be¬ 
tween the calculated positions from Eq. (6) and the undis¬ 
torted positions in ideal grid map. The two sets of coefficient 
a, and bj can be used for the corrections of all image signals. 

V. RESULTS AND DISCUSSIONS 

A. Effectiveness of the algorithm 

This algorithm was applied to the images acquired by our 
prototype lens-coupled digital x-ray imaging system. The 
corrected images in Fig. 6 showed a significant distortion 
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Fig. 6. Digital x-ray images showing the effect of distortion, (a) The cor¬ 
rected result of the grid image by fourth-order polynomials (see the image of 
the grid before the correction in Fig. 4). The mean residual error was re¬ 
duced from 8.04 pixel to 1.78 pixel, (b) A standard grid map (dashed lines) 
was overlaid on the corrected grid image to show the minimized geometrical 
distortion. Compare with Fig. 4 to see the effect of the correction. 


reduction when compared with the uncorrected images in 
Fig. 4. 

To quantitatively measure the effectiveness of the correc¬ 
tion algorithm, the absolute mean residual error (positional 
displacement of the grid image from its true location) was 
calculated for all of the intersections of the grid wires. The 
mean residual error before correction was 8.04 pixels and 
after corrections 1.78 pixels. In Fig. 7, the mean values of the 
radial distortion error, Ap, before and after correction are 
plotted as function of the field angle a of the optical system 
(refer to Fig. 1 for the definition of a). In Fig. 8, the mean 
values of the combined distortion error, Ay, before and after 
correction, are also plotted versus the field angle a (refer to 
Fig. 5 for illustration of Ap and Ay). 

Clearly, the polynomial algorithm effectively corrected 
the lens distortion in digital x-ray imaging. 



Fig. 7. The mean radial distortion errors before and after correction as a 
function of field angle of the optical system. Before the correction, the error 
ranges from 1.26 pixels (0.109 mm) to 7.42 pixels (0.64 mm). After correc¬ 
tion, the mean error ranges from 0.56 pixels (0.05 mm) to 1.00 pixels (0.086 
mm). 
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Fig. 8. The combined distortion errors Ay with respect to field angle before 
and after correction (refer to Fig. 5 for the definition of Ay). Before the 
correction, the error ranges from 2.39 pixels (0.206 mm) to 9.71 pixels (0.84 
mm). After correction, the error ranges from 1.45 pixels (0.125 mm) to 2.18 
pixels (0.19 mm). 

B. Order of polynomials and its impact on distortion 
correction 

As described previously, two sets of N-order polynomials 
were used in correcting both radial and tangential distortions. 
It is important to select a proper order of N to achieve opti¬ 
mized result. Theoretically, higher-order polynomial leads to 
less residual error. However, higher-order polynomials can 
also cause oscillation and make the mathematical solutions 
unstable. In our experiments, the effect of different orders of 
the polynomials on lens distortion correction was tested by 
generating sets of coefficients for orders: N= 2, 3, 4, 5, 6. 
The mean residual errors after correction for different N val¬ 
ues are listed in Table II. The results showed that the image 
correction improved as the order of the polynomial increased 
up to four. However, the correction effect decreased with 
higher-order polynomials. 

In addition to the possible oscillation, another reason for 
lack of improvement was that it could lead to some extra 
round-off errors when the higher order polynomial transfor¬ 
mation was performed. For the investigations using our im¬ 
aging system, the optimal order of polynomials appeared to 
be four. The corresponding polynomial coefficients a t and b t 
are listed in Table III. 

VI. CONCLUSION 

Optical coupling is an important link in various digital 
x-ray imaging systems. The optical properties of the compo¬ 
nents, either a lens, a fiber bundle or other optoelectronic 
device will affect the quality of clinical procedures signifi¬ 
cantly. There is no ideal optical component or optoelectronic 


Table II. The mean residual errors of corrected image using polynomial 
transformation of different orders. 


Polynomial order (N) 

2 

3 

4 

5 

6 

Mean residual error (Pixel) 

2.14 

2.02 

1.78 

3.94 

4.77 

Mean residual error (mm) 

0.18 

0.17 

0.15 

0.34 

0.41 


Table III. The polynomial coefficients used for correcting lens distortions 
(W=4). 


For radial 



*3 

"4 

Distortion 

1.008 

2.836X 1<T 5 

- 1.634X 1<T 7 

1.091X10’ 10 

For tangential 


*2 

*3 


Distortion 

1.003 

-2.293X10" 4 

-4.751 X10 -4 

3.074X 10" 5 


device that can relay an image without degrading its quality 
or distorting its shape. One could, however, balance the de¬ 
sign trade-offs of the selected components according to the 
specific clinical applications. 

Arguably, optical lens still is the most conventional and 
widely used imaging component. In comparison with many 
other optical techniques, lenses have many advantages, such 
as (1) cost effectiveness: It may be expensive to build just 
one or two custom lenses, but in mass production, lenses can 
be made at very low cost; (2) flexibility: Even with a simple 
relay lens, the parameters such as magnification, aperture, 
field of coverage, and depth of focus are adjustable; (3) im¬ 
age quality: A well-designed lens offers much higher spatial 
resolution than many other optical components or optoelec¬ 
tronic devices. For most commercial and custom lenses, ab¬ 
errations are selectively minimized according to the specific 
applications. However, residual aberrations, including geo¬ 
metrical distortions do exist even with a well-designed and 
well-built lens. Although these residual distortions are usu¬ 
ally small in terms of absolute or relative errors, they will 
reduce, for instance, the positional accuracy of image guided 
procedures. 

The distortion introduced by relay lenses, however, can be 
minimized with digital image processing techniques. In fact, 
it is easier to correct distortions introduced by lenses than 
that by many other types of imaging components, since lens 
distortion is primarily radial and tangential dependant. 24 27 
In other words, lens distortion may be “mapped” with a set 
of mathematical formulas and compensated accordingly us¬ 
ing appropriate algorithms, as demonstrated in this paper. 

We have shown that based on the nature and the classifi¬ 
cations of lens distortion, the coefficients of the polynomial 
in the correction algorithm can be determined to minimize 
both radial and tangential distortions. For the specific imag¬ 
ing system investigated, the fourth order polynomials pro¬ 
vided the best correction results. The low order polynomials 
were not sufficient for modeling the lens distortion, and the 
higher order polynomials tended to introduce oscillation that 
makes the algorithm unstable. 

Overall, our analysis suggests that lens distortion should 
be considered in developing optically coupled digital x-ray 
imaging systems for position-dependant, quantitative x-ray 
imaging, and that residual geometrical distortions could be 
minimized effectively using digital image processing tech¬ 
niques. The advantages of using lens to relay images from 
scintillating screen to electronic imager are cost effective¬ 
ness, flexibility and image quality. 


Medical Physics, Vol. 27, No. 5, May 2000 





















912 


Liu et a/.: Lens distortion 


912 


ACKNOWLEDGMENTS 

This work was supported in part by a NIH-NCI grant 
(Grant No. CA 70209) (P.I.: H Liu) and by a U.S. Army 
breast cancer research Grant No. DAMD 17-97-1-7138 (P.I.: 
H Liu). 

a) Author to whom correspondence should be addressed. Telephone: (410) 
614-9291; Fax: (410) 614-7760. Electronic mail: hliu@jhim.edu 
1 A. Karellas, L. Harris, H. Liu, M. Davis, and C. D’Orsi, “Charge- 
Coupled Device detector: Performance considerations and potential for 
mammographic imaging,’’ Med. Phys. 19, 1015-1023 (1992). 

2 H. Liu, J. Xu, G. Halama, and J. McAdoo, “Digital fluoroscopy using an 
optically coupled CCD: Preliminary investigations,’’ SPIE Proceedings 
2976, 256-261 (1997). 

3 A. Jalink, J. McAdoo, G. Halama, and H. Liu, “CCD-mosaic technique 
for large field digital mammography,’’ IEEE Trans. Med. Imaging 15, 
260-267 (1996). 

4 H. Liu, L. L. Fajardo, J. R. Barrett, and R. A. Baxter, “Contrast-detail 
detectability analysis: comparison of a digital spot mammography system 
and an analog screen-film mammography system,” Acad. Radiol. 4, 
197-203 (1997). 

5 H. Liu, L. L. Fajardo, and B. Penny, “Signal-to-noise ratio, noise power 
spectrum and detective quantum efficiency analysis of optically coupled 
CCD mammography imaging systems,” Acad. Radiol. 3, 799-805 
(1996). 

6 H. Liu, A. Karellas, S. Moore, and L. Harris, “Lesion detectability con¬ 
siderations for an optically coupled CCD x-ray imaging system,” IEEE 
Trans. Nucl. Sci. 41, 1506-1509 (1994). 

7 H. Liu, A. Karellas, L. Harris, and C. D’Orsi, “Optical properties of fiber 
tapers and their impact on the performance of a fiber optically coupled 
CCD x-ray imaging system,” SPIE Proceedings 1894, 136-147 (1993). 
8 H. Liu, A. Karellas, and L. Harris, “Methods to calculate lens coupling 
efficiency in an optically coupled CCD x-ray imaging system,” Med. 
Phys. 21, 1193-1195 (1994). 

9 B. H. Walker, “Primary lens aberrations,” in Optical engineering fun¬ 
damentals (SPIE Optical Engineering Press, Bellingham, Washington, 
1998), pp. 129-148. 

10 J. M. Boone, J. A. Seibert, W. A. Barrett, and E. A. Blood, “Analysis and 
correction of imperfections in the image-intensifier-tv-digitizer imaging 
chain,” Med. Phys. 18, 236-242 (1991). 
n R. Ning, J. K. Riek, and D. L. Conover, “An image intensifier based 
volume tomographic angiography imaging system: geometric distortion 
correction,” SPIE Proceedings 2163, 199-210 (1994). 


12 B. Achueler and X. Hu, “Correction of image intensifier distortion for 
three-dimensional x-ray angiography,” SPIE Proceedings 2432, 272-279 
(1996). 

13 E. Gronenschild, “The accuracy and reproducibility of a global method 
to correct for geometric image distortion in the x-ray imaging chain,” 
Med. Phys. 24, 1875-1888 (1997). 

14 E. L. Hall, Computer Image Processing and Recognition (Academic, 
New York, 1979), pp. 185-189. 

15 K. R. Castleman, Digital image Processing (Prentice Hall, Upper Saddle 
River, New Jersey, 1996), pp. 120-127. 

16 S. Shah and J. Aggarwal, “Intrinsic parameter calibration procedure for a 
(high distortion) fish-eye lens camera with distortion model and accuracy 
estimation,” Pattern Recogn. 29, 1775-1788 (1996). 

,7 Lasmea, “Camera Calibration from Spheres Images,” Proceedings of 
Computer Vision-ECCV’94, Stockholm, Sweden (Springer-Verlag, Ber¬ 
lin, 1994), pp. 449-454. 

18 D. C. Brown, “Close-range camera calibration,” Photogramm. Eng. Re¬ 
mote Sens. 37, 855-866 (1971). 

19 Manual of Photogrammetry, 4th ed. (Amer. Soc. Photogrammetry, Falls 
Church, Virginia, 1980). 

20 D. Ballard and C. Brown, Computer Vision (Prentice-Hall, Inc., Engle¬ 
wood Cliffs, New Jersey, 1982). 

21 J. Weng, P. Cohen, and M. Hemion, “Camera calibration with distortion 
models and accuracy evaluation,” IEEE Trans. On PAMI-14, No. 10, 
965-980 (1992). 

22 Z. Jericevic, D. M. Benson, J. Bryan, and L. C. Smith, “Geometric cor¬ 
rection of digital images using orthonormal decomposition,” J. Microsc. 
149, Pt.3 233-245 (1988). 

23 R. Fahrig, M. Moreau, and D. M. Holdsworth, “Three-dimensional com¬ 
puted tomographic reconstruction using a C-arm mount XRII: correction 
of image intensifier distortion,” Med. Phys. 24, 1097-1106 (1997). 

24 D. C. Brown, “Decentering distortion of lenses,” Photogramm. Eng. 
Remote Sens. 32, 444-462 (1966). 

25 R. Y. Tasi, “A versatile camera calibration technique for high-accuracy 
3D machine vision metrology using off-the-shelf TV cameras and 
lenses,” IEEE Journal of Robotics and Automation 3, 323-344 (1987). 

26 Y. Nomura, M. Sagara, H. Naruse, and A. Ide, “Simple calibration algo¬ 
rithm for high-distortion-lens camera,” IEEE Trans. On PAMI-14, No. 
11, 1095-1099 (1992). 

27 H. Bacakoglu and M. S. Kamel, “A three-step camera calibration 
method,” IEEE Trans. Instrum, and Measurement 46, 1165-1172 
(1997). 


Medical Physics, Vol. 27, No. 5, May 2000 






Noise power spectra of images from digital mammography detectors 
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Noise characterization through estimation of the noise power spectrum (NPS) is a central compo¬ 
nent of the evaluation of digital x-ray systems. We begin with a brief review of the fundamentals of 
NPS theory and measurement, derive explicit expressions for calculation of the one- and two- 
dimensional (ID and 2D) NPS, and discuss some of the considerations and tradeoffs when these 
concepts are applied to digital systems. Measurements of the NPS of two detectors for digital 
mammography are presented to illustrate some of the implications of the choices available. For both 
systems, two-dimensional noise power spectra obtained over a range of input fluence exhibit pro¬ 
nounced asymmetry between the orthogonal frequency dimensions. The 2D spectra of both systems 
also demonstrate dominant structures both on and off the primary frequency axes indicative of 
periodic noise components. Although the two systems share many common noise characteristics, 
there are significant differences, including markedly different dark-noise magnitudes, differences in 
NPS shape as a function of both spatial frequency and exposure, and differences in the natures of 
the residual fixed pattern noise following flat fielding corrections. For low x-ray exposures, quan¬ 
tum noise-limited operation may be possible only at low spatial frequency. Depending on the 
method of obtaining the ID NPS (i.e., synthetic slit scanning or slice extraction from the 2D NPS), 
on-axis periodic structures can be misleadingly smoothed or missed entirely. Our measurements 
indicate that for these systems, 1D spectra useful for the purpose of detective quantum efficiency 
calculation may be obtained from thin cuts through the central portion of the calculated 2D NPS. 

On the other hand, low-frequency spectral values do not converge to an asymptotic value with 
increasing slit length when ID spectra are generated using the scanned synthetic slit method. 
Aliasing can contribute significantly to the digital NPS, especially near the Nyquist frequency. 
Calculation of the theoretical presampling NPS and explicit inclusion of aliased noise power shows 
good agreement with measured values. © 1999 American Association of Physicists in Medicine. 
[S0094-2405(99)00707-5] 

Key words: noise power spectrum, digital mammography, detector, image analysis 


I. INTRODUCTION 

The signal-to-noise ratio imposes the fundamental limitation 
to object perceptibility in a digital radiograph because image 
contrast can be manipulated during the display of digitally 
acquired radiographic images. Thus, noise characterization 
plays an increasingly central role in the evaluation of system 
performance in medical imaging. Well-developed methods 
used for noise characterization in film-based systems are cur¬ 
rently being reassessed and reformulated to accommodate 
both the different nature of the noise and to take advantage 
of the analytical tools available with digital processing. An 
indication of both the current interest level and lack of con¬ 
sensus regarding noise characterization is the recent estab¬ 
lishment of an AAPM Task Group (No. 16) for Standards for 
Noise Power Spectrum Analysis. 

The noise power spectrum (NPS) is a spectral decompo¬ 
sition of the variance. As such, the NPS of a digital radio- 
graphic image provides an estimate of the spatial frequency 
dependence of the pixel-to-pixel fluctuations present in the 
image. Such fluctuations are due to the shot (quantum) noise 
in the x-ray quanta incident on the detector, and any noise 
introduced by the series of conversions and transmissions of 
quanta in the cascaded stages between detector input and 
output. Examples of the latter are gain variances in the con¬ 


version of x-ray quanta to light quanta in a phosphor (or to 
electron-hole pairs in a solid-state detector), statistical fluc¬ 
tuations in the transmission of optical quanta between a scin¬ 
tillator and a photodetector, and additive noise sources such 
as preamplifier noise. The NPS is a much more complete 
description of image noise than is quantification of integrated 
(total) noise via simple measurement of the rms pixel fluc¬ 
tuations, because it gives information on the distribution in 
frequency space of the noise power. An understanding of the 
frequency content of image noise can provide insight regard¬ 
ing its clinical impact. For example, in mammography, ex¬ 
cess high-frequency noise may render the detection of micro¬ 
calcifications impossible. If desired, the total variance can be 
obtained by integrating the NPS over spatial frequency. 

In this paper we discuss aspects of NPS estimation as 
applied to digital radiographic systems. We begin with a 
general review of pertinent definitions and relations, and dis¬ 
cuss the most important considerations when applying spec¬ 
tral estimation techniques to digital radiographic systems. 
Mathematical derivations sufficient to provide the reader 
with tools for quantitative application of the concepts are 
given, with more detailed development available in the cited 
references. These techniques are then applied to the analysis 
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of the noise properties of two developmental detectors for 
full-field digital mammography using charge coupled de¬ 
vices (CCDs). 

II. BACKGROUND 

Linear systems analysis can be usefully applied to any 
radiographic system whose transfer characteristic between 
input and output is linear (or linearizable) and stationary (the 
transfer properties are independent of location at the input 
surface). 1,2 In this formalism, signal and noise are typically 
decomposed into their Fourier components (although other 
bases can be used), and input and output relations are ex¬ 
pressed in terms of the frequency space Fourier amplitudes. 
Digital radiographic systems such as the one considered here 
usually exhibit linear responses to the incident x-ray fluence 
over several orders of magnitude. The following analysis as¬ 
sumes that the detector is operated within its linear range. 
For the digital mammography systems whose noise charac¬ 
terization is reported here, care was taken to first determine 
the linear range of operation for each. The range of expo¬ 
sures used for the NPS estimates in this paper fall well 
within the boundaries of linear operation. In a digital detec¬ 
tor the analog (presampling) signal and noise propagated 
through the detector are sampled by the digital matrix at 
discrete points. Therefore, the system response is not strictly 
shift invariant (unless shifts are by an integer number of 
pixels). The effects of digital sampling on the NPS are de¬ 
scribed below. 

A random process such as the noise fluctuations in an 
image is considered (wide-sense) statistically stationary , if 
its autocorrelation function (see below) is independent of the 
particular data sample (e.g., the particular location in the 
image) used to obtain it. This is only strictly true in a radio- 
graphic image if (a) noise samples are obtained from spa¬ 
tially uniform exposures, and (b) the detector adds no spa¬ 
tially fixed correlated noise. Dobbins has shown that the 
presence of a deterministic signal in the image data used to 
estimate the NPS has no effect on the resultant spectral 
estimates. 3 One other important assumption concerning the 
noise is whether or not it is ergodic. The noise in an image is 
ergodic if expectation values (averages) obtained from data 
samples at various locations in the image are equivalent to 
ensemble averages obtained from repeated measurements un¬ 
der identical conditions at a single location. Ergodicity im¬ 
plies stationarity, but not necessarily vice versa. The assump¬ 
tion of ergodicity of the noise in radiographic images made 
under conditions of uniform irradiation is typically made, 
and permits averaging of spectra obtained from several areas 
of the image. Strictly speaking, even under conditions of 

_I 


nominally uniform irradiation, factors such as the heel effect 
result in deviation from true ergodicity. For this reason, it is 
often best to estimate the NPS from a common, limited area 
of a series of images, rather than from a large portion of a 
single image. In the work described here, data used for a 
given spectral estimate were obtained from a series of im¬ 
ages, but using an area located at a fixed distance from the 
chest wall edge. 

In a cascaded system such as a digital mammography de¬ 
tector, the output of one stage constitutes the input to the 
following stage. In a linear, shift invariant cascaded system, 
each of the independent stages are also linear and shift in¬ 
variant. At the output of any given stage of a cascaded de¬ 
tector system, the minimum possible value of the NPS is that 
set by the (uncorrelated) Poisson variance of the output 
quanta, which is equal to the number of quanta. Building on 
the work of Rossmann 4 and Kemperman and Trabka, 5 Rab- 
bani, Shaw, and Van Metter have developed expressions for 
the propagation of the NPS through cascaded linear 
systems. 6 Using moment generating functions, they show 
that noise generated during a stochastic amplification process 
in one stage (such as the conversion of a single absorbed 
x-ray photon to many light photons) is effectively passed as 
signal to the subsequent stage. Furthermore, they show that 
in a stochastic scattering process, while correlated compo¬ 
nents of the NPS are reduced by the square of the modulation 
transfer function (MTF) associated with the blurring process, 
uncorrelated Poisson noise is transferred independently of 
the MTF. This stochastic type of blur has a different effect 
on the NPS than deterministic blur, in which both correlated 
and uncorrelated components of the NPS of the preceding 
stage are attenuated by the square of the blur MTF. 7 Deter¬ 
ministic blur does not degrade the signal-to-noise ratio be¬ 
cause the signal is also modulated by the blur MTF. Visible 
photon scatter in a phosphor is an example of stochastic blur; 
the effect of the integrating aperture of the individual pixels 
in a digital system is an example of deterministic blur. 

A. Direct and indirect methods of NPS calculation 

The NPS is defined as the Fourier transform (FT) of the 
autocorrelation function, defined in one dimension as 8 

1 CL/2 

C{x)= lim- p{x+T)p*{T)dT, (1) 

I_00 L j -LI2 

where p(x) is the value (pixel value for a digital image) of 
the one-dimensional image (ID) at position x , and p*(x) is 
its complex conjugate. Since p(x) is real,/?( jc) =/?*(*). The 
Fourier transform of C(x) is then 


foe 1 fU2 foe 

S(u)= C(x)e~ 2mxu dx = lim — /?*(r)e + 27n ™ p(x+T)e~ 2mxu dx[e~ 2iriru ]dT 

J -oo L-+ oc^ J-L/2 J-oo 


1 f L/2 1 fL/2 

= lim— p*(T)e^ 2niTU p(a)e~ 2ni<ru dadr= lim— p*(T)e +2mTU dTP(u). (2) 

J-L/2 J -oo 
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Inspection of Eqs. (1) and (2) shows that 

S( M )-FT{C(x)}=limj|/>(u)| 2 . (3) 

Thus, the NPS may be calculated either from the Fourier 
transform of the autocorrelation function (the indirect 
method), or from the square of the modulus of the Fourier 
transform of the data itself (the direct method). Note that 
since the data p(x) are real, that P(-u) = P*(u). That is, 
the NPS at a given negative frequency is equal in magnitude 
to the NPS at the corresponding positive frequency. Extend¬ 
ing the expression of Eq. (3) to two dimensions: 

$(«,») = FT{C(*o>)}= lim r^\P(u,v)\ 2 . (4) 

A\y->o o AI 

With the advent of the fast Fourier transform (FFT) and fast 
computers, indirect calculation of the NPS via the autocorre¬ 
lation function has largely been replaced by the direct 
method. This paper will discuss only the latter method. 

If the stationary random process being characterized by 
the NPS is also ergodic, as is typically the case for radio- 
graphic image noise, then the spectral estimate is found by 
ensemble averaging. That is, the final NPS is the average of 
spectra obtained from a series of uniform irradiation images. 
Then, Eq. (4) becomes 

I i I ext 2 cm 

S(u,v)= lim — p(x,y) 

X,Y-k*\ AY I J -X/2J -Y/2 

Xe -2*i(ux+vy) dxdy ^ ^ ( 5) 

where the pointed brackets denote ensemble averaging. 


B. One-dimensional NPS 


The detective quantum efficiency (DQE) is perhaps the 
best single descriptor of radiographic detector performance. 
It quantifies how well the detector is able to transfer the 
signal-to-noise ratio (SNR) inherent in the x-ray fluence at its 
input to the resulting image. For most two-dimensional (2D) 
digital imaging detectors, the DQE is a function of at least 
three independent variables: its two spatial coordinates, and 
the x-ray exposure. It is a function of a single spatial coor¬ 
dinate for the special case of detectors with an isoplanatic 
response. The DQE can be written: 


DQE (w,u,<p) = 


(<pg) 2 C(M,i;)MTF 2 (u,u) 

■S(w,u,<p)SNR 2 n 


(6) 


where <p represents the magnitude of the input x-ray fluence 
at the detector surface; g relates changes in the zero- 
frequency output signal to that at the input (i.e., it is the 
system gain); C(u,v) is the Fourier decomposition of the 
input signal, normalized to unity at maximum amplitude; 
MTF(w,u) is the detector modulation transfer function; and 
S(w,i>,<p) is the output noise power spectrum. 9 Although 
5(w,u,<p) is a function of the input fluence <p, for notational 
simplicity it will be referred to hereafter simply as S(u,v). 


For a monochromatic input x-ray spectrum, or in the case of 
a photon-counting x-ray detector, SNR in is given by the mag¬ 
nitude of the input fluence <p. 

In part because illustration of DQE (u,v,<p) requires a 
four-dimensional medium, DQEs of 2D imaging systems are 
typically presented one spatial frequency coordinate at a time 
for clarity. By definition, 

MTF(w,u) = |2D FT{psf(x,j')}|, (7) 

where psf(jc,_y) is the spatial response of the detector to a 
point x-ray input, or point spread function, and 2D FT rep¬ 
resents a two-dimensional Fourier transform. The ID MTF is 
typically measured using a narrow slit or edge, 10 and is 

MTF(w) = FT{lsf(x)} = FT j ^ [esf(x)] J, (8) 

where lsf(x) = /“ocPsf {x,y)dy is the line spread function and 
esf(x) is the edge spread function. Then, 

J \sf(x)e- 2nixu dx=\ I * {ps f(x,y)dy}e- 2irixu dx 

J -00 j J — oc 

= |2D FT{psf(x,>0}Uo 

= MTF(w,0). (9) 

Equations (8) and (9) indicate that evaluation of FT{lsf(^)} 
yields MTF(w,0), and thus the appropriate expression of the 
NPS for calculation of the single coordinate DQE is S(ufi). 
Similar arguments pertain to MTF(y). n 

Up to now, most characterizations of S(w,0) [or 5(0,i;)] 
of digital radiographic systems have been made using the 
synthesized slit method. 12-15 In the original development of 
the method in the analog (film-based) context, the optical 
aperture of the microdensitometer used to measure the film 
optical density formed the slit. The methodology has since 
been applied to digital systems by forming an effective slit 
from rectangular groups of contiguous pixels. In either case, 
a long, narrow slit is scanned in steps in one direction over 
the image. At each step the average intensity in the slit (av¬ 
erage transmitted light for a scanning microdensitometer, or 
average pixel value for a digital system), is recorded. In this 
way, a one-dimensional data series is generated. The squared 
modulus of the Fourier transform of the series, after scaling, 
gives the measured ID NPS. In the limit of an infinitely long 
ID data series in the x direction, the resulting NPS is 16,17 

Sss(“) = /_ (10) 

where T(u,v) is the optical transfer function (OTF) of the 
slit, and where S ss (u) is the ID NPS measured via the 
scanned slit method. If the slit is rectangular, with 
x-dimension w and y-dimension Z,, then its point spread func¬ 
tion is 

A ^ )= (i rect ^)(r ect z)’ (11) 

and therefore, T(u,v) is 
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r(w,u) = FT{/i(x,y)} = sinc(wM)sinc(Z,i;), (12) 

where the sine function is defined as 
sin 7 tjc 

sinc(x) =-. (13) 

7 TX 

Equation (10) becomes 

S ss (u) = sinc 2 (wu) J S(u,v)sinc 2 (Lv)dv. (14) 

if the slit width is chosen sufficiently narrow, then 
sinc(ww)~l within the range of u over which S ss (u) is to 
be evaluated. In any case, correction of the measured NPS by 
this factor is straightforward. In a digital system, a synthe¬ 
sized slit of 1 X N l pixels constitutes sampling the presam¬ 
pling NPS (the NPS smoothed by the nonzero pixel aperture, 
but prior to sampling at discrete intervals given by the pixel 
spacing), using a line of N L delta functions, at an interval in 
the scan direction equal to the center-to-center pixel spacing, 
Ax. The N l sampled values are then averaged. Thus, 
sinc(ww)—>sinc(0) = 1, and L = N L Ay, where A y is the 
pixel size perpendicular to the scan direction. 

As L—► oo, sinc 2 (Z,u) becomes sufficiently narrow that 
S(u,v) is approximately constant over the range of v (near 
v = 0) in which sinc 2 (Iu) is appreciably nonzero, and 

S ss ( W ) = J S(u,v)s\nc 2 (Lv)dv 

f* S(u y 0) 

=S(w,0) sinc 2 (Z,u)^u=—-—. (15) 

J -00 L 

Thus, the desired NPS is S(w,0) = Z,S ss (w). 

Substantial effort has been devoted to understanding the 
criteria for “sufficiently large L.” Sandrik and Wagner stud¬ 
ied the effect of varying slit length on the 1D NPS of radio- 
graphic screen-film systems. 18 They found that the magni¬ 
tude of the low-frequency NPS is underestimated if the slit 
length is insufficient, but that the low-frequency components 
approach a limiting value as the slit length is increased. The 
slit length required for low-frequency NPS values within 5% 
of the plateau value was dependent on the particular screen- 
film combination, but ranged between a value equal to the 
length of the scan and one nearly twice as long. Koedooder, 
Strackee, and Venema developed an expression to numeri¬ 
cally approximate the integral in Eq. (14). ig The technique is 
applicable only to systems with rotational symmetry. 

An alternate to the scanned slit approach is to extract a 
slice along one of the primary spatial frequency axes from 
the 2D NPS. However, because these slices contain noise 
power from the orthogonal dimension near to and at zero 
frequency, low-frequency trends in the data can contribute a 
significant portion of the noise power unless these trends are 
first removed. For this reason, the limited number of studies 
employing this approach have used “thick” slices made by 
extracting NPS data near to, but not directly on, the primary 
axes. 20 21 As will be demonstrated below, in some digital 
detectors this approach can overlook structures that occur 
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only at zero frequency in the orthogonal dimension (i.e., they 
are constant in the orthogonal spatial direction, such as a 
stripe or band). 


C. Digital noise power spectrum 

In the above expressions, the 1D and 2D NPS are written 
as functions of a continuous variable p(x) that describes the 
point-to-point fluctuations in the image. In a digital system, 
this continuous function of the spatial coordinates is sampled 
at regular intervals, Ax. Thus, p{x) is evaluated at a set of 
discrete locations, x = n Ax, n = 0,1,2, ...JJ, where X=NAx. 
Correspondingly, functions of the spatial frequency variable 
u are now evaluated at discrete frequencies given by u 
= &Aw, k=0,± 1,±2,.... The maximum spatial frequency 
sampled, the Nyquist frequency, is u N = 1/(2Ax). To experi¬ 
mentally estimate the digital NPS, the above continuous ex¬ 
pressions must be rewritten in terms of discrete variables. 
With the preceding definitions, Eq. (3) can be written 

S d (u) = S(kAu) 

^ N-\ 2 

= Iim ttt- Ax 2 p(n\x)e- 2lriU, * u)(n£ix) 


= lim 

N&X—>ao 


Ax 

nr 


N- 1 

2 p(n&x)e 

n = 0 


2 


2iri(klu)(nAx) 


(16) 


where N is the number of samples in the interval X and the 
variable u now takes on the discrete values 0, =ti 1/A^, 
± 2/X ,.... Similarly, Eq. (4) for the 2D NPS becomes 


S d (u 9 v) = S(jAu 9 kAv) 


= lim 

N&x —*« 

M&y—tcc 


Ax Ay 

nr ~m 


N-\ M -1 

X 2 p(«Ax,mA^) 

n = 0 m = 0 


2 


X e -2Tri[(Jbu)(nbx) + (kbv)(m&y)] 


(17) 


The experimentally determined estimate of the digital NPS 
differs from the underlying analog spectrum in the following 
ways: 

(a) The noise has been convolved with the rect (boxcar) 
function given by the pixel aperture (deterministic blurring). 

(b) The noise data has been sampled at discrete intervals 
given by the pixel center-to-center spacing (equivalent to 
multiplication by a Comb function). 

(c) The record (image data) used to compute the estimate 
has been truncated from its original infinite length to a length 
given in each dimension by the number of points times the 
sampling interval (equivalent to multiplication by a 2D box¬ 
car function of dimensions A/Ax X A/Ay). 

The result of (a) is the multiplication of the analog NPS 
with the absolute square of the 2D sine function that is the 
Fourier transform of the pixel aperture function. The result is 
often referred to as the presampling NPS: 

Spre( w >t;) = S(w,i;)|<z6 sinc(aw)sinc(Au)|, (18) 
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where S(u,v) is the analog NPS, and a and b are the pixel 
aperture dimensions. 

The result of (b) is that the presampling NPS is convolved 
with the transform of the sampling Comb function, 
III(w,i;;(l/Ax)(l/Ay)): 


S d («,u) = Sp re (M,u)®®IIl|w,w;^^j, (19) 

where Ax and A y are the sampling intervals in the x and y 
dimensions, and <8>® denotes two-dimensional convolution. 
Note that this Comb function consists of an infinite array of 
delta functions separated in each dimension in frequency 
space by two times the respective Nyquist frequencies; U N 
= l/(2Ax), u N =l/(2Ay). S d (u,v) is, therefore, an infinite 
2D array of replications of S pre (w,i;), with the origin of each 
S pre separated by 2 u N in the u direction and 2v^ in the v 
direction. This replication means that noise beyond the Ny¬ 
quist frequency overlaps noise below the Nyquist frequency 
of adjacent replications. This overlap, or aliasing , becomes 
increasingly severe as the spatial sampling interval is in¬ 
creased, that is, as S pre is increasingly undersampled. The 
presampling NPS cannot be directly measured using fine 
sampling techniques such as those commonly employed to 
measure the presampling MTF, because the phases of the 
Fourier components of the image noise are random. 22 

Using digitized screen-film images, Giger, Doi, and Metz 
have studied the effects on the digital NPS of the sizes of the 
integrating pixel aperture and the sampling interval. 23 Their 
results show that aliasing of high-frequency noise (due, for 
example, to the high-frequency components of x-ray quan¬ 
tum noise) above the Nyquist frequency can form a signifi¬ 
cant component of the measured digital NPS. Large integrat¬ 
ing aperture size (i.e., larger than the sampling interval) 
tends to reduce the effect of this aliasing by smoothing the 
high-frequency x-ray quantum noise, at the expense of spa¬ 
tial resolution. Aperture sizes less than the sampling distance 
increase the overall noise level, and particularly so at high 
frequencies. Dobbins has recently discussed the impact of 
undersampling on the interpretation of the MTF of digital 
systems. 3 

The result of (c) is that the NPS is also convolved with the 
square of the transform of the data window w(x ,y): 

S dtW (u,v)=S/,u,v)®®\W(u,v)\ 2 , (20) 

where S d (u,v) is the measured NPS, and the frequency win¬ 
dow, W(u,v) is the Fourier transform of the data window 
w(x,y). 24 

If the data series is simply truncated after NXM points, 
corresponding to multiplication of the infinite data sequence 
by a unity height box of dimensions VAxXA/Ay, then the 
NPS is convolved with a frequency window with modulus 


\W(u,v)\ = NbxMby 


sin( ttN&x u ) sin( irM A yv) 
Tr 2 NAxuM Ayv 


( 21 ) 


A variety of data windows have been studied in many 
signal processing contexts, with the optimum choice for a 


given situation depending on the characteristics of the par¬ 
ticular spectrum being characterized. Selection of a data win¬ 
dow is based on the shape and width of the resulting fre¬ 
quency window. In general, there is a compromise between 
the width of the central lobe and the magnitude and extent of 
any side lobes. A broad central lobe reduces spectral resolu¬ 
tion, and significant side lobes produce leakage , wherein 
noise power from nearby frequencies are mixed in the con¬ 
volution process. In the case of power spectra exhibiting 
strong periodic components, as is the case in many digital 
radiographic systems, minimization of side lobe magnitude 
is often necessary for distinguishing small amplitude features 
adjacent to large amplitude peaks. 25,26 


D. Experimental estimation of the noise power 
spectrum 


1. Uncertainty in the spectral estimates 


Because noise in a radiographic image arises from sto¬ 
chastic (random) processes, any realization (image) of the 
noise associated with those processes is merely a sample 
chosen from the ensemble of all possible samples. Thus, 
there is an inherent statistical uncertainty in the representa¬ 
tion of the underlying (true) population variance by that 
sample. In this regard, Fourier decomposition of image noise 
(such as is performed in determination of the NPS) is differ¬ 
ent from the decomposition of signal (such as is performed 
in determination of the modulation transfer function), where 
amplitudes and phases are known and fixed. The variance in 
the spectral components of the measured NPS is independent 
of the number of data points (e.g., the number of pixels) in 
the sample used to calculate it, and for mean-zero, normally 
distributed noise, results in a coefficient of variation (COV) 
of approximately unity for each spectral estimate, S(u,v) 
except at zero frequency and the Nyquist frequency, where 
the COV is 2 1/2 . 27 The variance (uncertainty) in the estimate 
can be reduced by sectioning the entire VAxXA/Ay noise 
record into a number of contiguous blocks, and averaging the 
spectra obtained from each (ensemble averaging). The COV 
of the final estimate decreases as the inverse square root of 
the number of determinations. However, the spectral resolu¬ 
tion is inversely proportional to the number of data points 
used to generate each NPS. Thus, for a fixed amount of data, 
there is a tradeoff between frequency resolution and variance 
in the final NPS estimate. For example, if N data points, 
spaced at an interval Ax, are used to calculate a ID NPS, 
then the spectral resolution and COV are, respectively, 27 


a/= a^’ C0Vs1 °- 


( 22 ) 


If these N points are divided into n successive sections, then 
the frequency resolution and COV become 


A/= 777\ = 777 —, CO V= -=. 

(£)*, NAx ^ 
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Table I. Summary of the various components of the NPS and methods for isolating each. The NPS components 
are (a) random (stochastic), (b) nonstochastic but varying from one image to another, and (c) nonstochastic and 
fixed from image to image (fixed pattern noise). 


NPS 

Method 

Components present 

sjtw) 

Average spectra from individual 
images 

Stochastic, varying nonstochastic, 
fixed pattern 


Average 0.5 times spectra from 
difference images 

Stochastic, varying nonstochastic 

S„ g (u,v) 

Spectrum from average of many 
images 

Varying nonstochastic, fixed pattern 

Sf. p .(u,v) 


Fixed pattern 

Sy,. n.(W,y) 

Savg( W , U ) — S f p ( U,V) 

Varying nonstochastic 

SstochCw.u) 

Su>t(u>v)-S„ g (u t v) 

Stochastic 


2. Stochastic and nonstochastic components 

Image noise can be categorized as (a) random (stochas¬ 
tic), (b) nonstochastic but varying from one image to an¬ 
other, or (c) nonstochastic and fixed from image to image 
(fixed pattern noise). For a complete description of image 
noise, each component should be analyzed. 

Fixed pattern noise usually arises from spatial variations 
in detector sensitivity. Although theoretically fixed pattern 
noise should be removed during routine postacquisition im¬ 
age processing (flat fielding), the removal may be incomplete 
or its completeness may be signal dependent. Note that since 
flat fielding typically involves division of the raw image by a 
normalized, highly averaged raw flood image, random noise 
in areas of low gain (average output signal per incident x 
ray) is amplified in the flat fielding process, while that in 
areas of high gain is reduced. Nonstochastic noise that is 
present intermittently or varies in location from image to 
image will not be reliably removed during flat fielding. An 
example of this type of noise is 60 Hz power line pickup 
during detector readout. 

The NPS obtained by averaging the spectra of many 
single uniform irradiation images, S tot (u,v ), contains compo¬ 
nents from all three types of noise. The NPS using the dif¬ 
ference of two images obtained under identical conditions of 
uniform irradiation (and divided by 2), 5 diff (w,t;), contains 
variance from stochastic and varying nonstochastic compo¬ 
nents, but none from fixed pattern noise. The magnitude of 
the noise power from fixed pattern noise can be determined 
from Sf p (u > v) = S tot (u,v)—S di f^u y v). Note that the NPS ob¬ 
tained from a highly averaged image, 5 avg (w,i;), differs from 
Sf.p.(u,v) by an amount equal to the noise power due to any 
time-varying nonstochastic components. The stochastic NPS 
components can be found by subtracting 5 avg (w,u) from 
5 tot (w^). 28 Table I summarizes the methods for isolating 
spectral components due to these three types of noise. 

III. METHODS 
A. Data acquisition 

Both mammographic systems evaluated utilize molybde¬ 
num target x-ray tubes, with molybdenum or rhodium filtra¬ 


tion. Data were acquired using Mo/Mo spectra and a variety 
of kVp and mAs settings. A block of acrylic 3.8 cm thick 
and sufficiently large to cover the entire detector surface was 
placed on the breast support to create realistic spectral com¬ 
position and scatter conditions. In selecting areas of the im¬ 
age for estimating noise power spectra, care was taken to 
avoid regions near the image periphery, where the scatter 
contribution is nonuniform. Detector entrance exposures 
were calculated by raising the acrylic block slightly above 
the detector surface, and using a thin ion chamber to measure 
the exposure at its exit surface. Corrections were then ap¬ 
plied to account for the slight difference between the location 
of the chamber and the detector entrance surface, relative to 
the focal spot. For each selected kVp and mAs, eight images 
were obtained. Data used for spectral estimates were taken 
from —5cmX5cm regions, located centrally in the image. 
The data regions were divided into the maximum number of 
contiguous sections consistent with the NPS estimation tech¬ 
nique to be used (e.g., NXL pixel sections for a scan in the 
x direction, or NX M pixels for a 2D NPS calculation). Spec¬ 
tra were calculated from each section of each of the eight 
images, and were then averaged. The final ensemble- 
averaged NPS estimates were thus averages of —300-1200 
spectra, resulting in standard errors of — 3%-6%. 

B. Subtraction of background trends 

Low-frequency background trends, such as those from the 
x-ray source heel effect, are often present in the data. Al¬ 
though such trends may have significant Fourier components 
only at frequencies below the lowest measurable frequency 
in the power spectrum (l IN Ax) 9 leakage into low-frequency 
spectral components can artificially inflate the low-frequency 
NPS. (As discussed in Sec. IIC, leakage is the spreading of 
spectral power from one frequency to adjacent frequencies, 
and is a result of the nonzero width of the frequency win¬ 
dow, which is the Fourier transform of the finite data win¬ 
dow.) For this reason, noise power spectra are frequently 
obtained by first subtracting one uniform exposure image 
from another made with the same exposure, and dividing the 
NPS, resulting from the subtracted image, by 2. The factor of 
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2 accounts for the fact that the variance in the subtracted 
image gets approximately equal contributions from the un¬ 
correlated noise in the two original images. Alternately, a 
difference image may be obtained by subtracting a highly 
averaged uniform exposure image from a single image. The 
random noise in the highly averaged image is negligible if a 
large number of images is used to obtain it. 

The image subtraction technique has the drawback, how¬ 
ever, of also removing any uncorrected fixed pattern noise 
occurring at higher frequencies. Thus, this method may not 
be optimum for estimating the NPS of clinical images fol¬ 
lowing uniformity correction (flat fielding), because noise 
removed in the subtraction process would nevertheless be 
present in clinical images. We have, therefore, evaluated two 
methods for removing low-frequency background trends 
from single images: subtraction of a low-pass-filtered version 
of the image, and subtraction of a low-order polynomial fit to 
the image. When using the first method, for simplicity, a 
rectangular (boxcar) filter was used to filter the entire 2D 
data set. Subtraction of the resulting filtered image from the 
original is equivalent to subtraction of a running local means 
as has been employed for one-dimensional data by others. 23 
Convolution in coordinate space with a rectangular filter is 
equivalent to multiplication in frequency space by a sine 
function. The dimensions of the square filtration kernel were 
chosen sufficiently large that the first zero crossing of the 
resulting sine function was at or below the minimum fre¬ 
quency sampled in the NPS. In the second method of back¬ 
ground trend subtraction, a two-dimensional polynomial was 
first fit to each uniform exposure image, and subtracted. Both 
first- and second-degree polynomials were evaluated. 

In order to compare the above two methods of back¬ 
ground trend subtraction, a synthetic data set was created. 
The data consist of a dc (average) pixel value, a normally 
distributed random component, a periodic component (one in 
each dimension), a linear (ramp) component, and a second- 
order component. The digital value p(i,j) of the pixel in the 
ith column and yth row was computed as 

p(ij)=d*[taad(i,j)]+d 2 (ij), (24) 

where the rand() operator generates one of a set of normally 
distributed random numbers with average zero and standard 
deviation one, d\ is a scaling factor, and 

d 2 (iJ) = a 0 sin(iu p +S u ) + b 0 sin(Jv p +S v ) 

+ a i i + b J + a 2 i 2 + bj 1 + C 0 , (25) 

where u p and v p are the spatial frequencies of the periodic 
components in the x and y directions, respectively, and S u 
and S v are phase factors. The constant C 0 sets the magnitude 
of the dc component, and the constants a 0 , b 0 , a x , b \, a 2 , 
and b 2 set the magnitudes of the periodic, linear, and qua¬ 
dratic components in the x and y directions, respectively. 
Figure 1 is a gray-scale plot of a portion of one of the syn¬ 
thetic data sets. 

Figure 2 shows a semilog plot of four representative ID 
power spectra generated from synthetic data. The spectrum 
labeled with diamonds is that of a data set with no back- 



Fig. 1. Gray-scale image of a synthetic data set used to evaluate methods of 
background trend removal. For ease of visualization, the data set shown 
contains only sinusoidal, random, and dc components (no linear or quadratic 
component). 


ground trend (a { = a 2 = b { = b 2 = 0). The second spectrum, 
labeled with squares, is from the same data, but with linear 
and quadratic components added. The third and fourth spec¬ 
tra (triangles and crosses), are calculated following trend re¬ 
moval via subtraction of a least-squares surface fit to a 
second-order polynomial or by subtraction of a boxcar- 
filtered version of the image, respectively. In this example, 
the boxcar filter was 64 X 64 pixels, the same size as the 
sections of data used to calculate the 2D NPS. The spectrum 
obtained with no trend removal demonstrates substantial 
leakage of low-frequency noise power into the frequency 
range of the periodic component. Figure 2 shows that both 
trend removal techniques reduce the low-frequency noise 
power sufficiently that it is no longer apparent in the spec¬ 
trum, even though the leakage is unchanged. The quadratic 
fit technique virtually restores the original NPS, a result not 
surprising given the quadratic nature of the trend. 

Figure 3 also shows the effects of different trend removal 
techniques, this time using flat-field image data obtained 
from one of the mammography detectors. The 2D spectra 
were computed from the average of 648 128X 128 pixel sec¬ 
tions. Thin slices were extracted along the v axis. Again, 
slices are shown from 2D spectra following no trend re¬ 
moval, and following subtraction of either a fitted quadratic 
surface or a low-pass-filtered data set. In addition, two dif¬ 
ferent size convolution kernels for low-pass filtering are 
compared. For a rectangular, nXm pixel convolution kernel 
of amplitude (wAx/wAy) -1 , where Ax and Ay are the pixel 
size in the x and y dimensions, the trend-corrected NPS is 
given by 

Sd t w,b = Sd,w[ 1 — sinc(/i Axw)sinc(wAyi;)]. (26) 

In order that only noise power close to dc is removed, one 
would like the term sinc(rtAxw)sinc(mAyi;) to be appre¬ 
ciable only for very low frequency. Therefore, the quantities 
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Fig. 2. Comparison of background trend removal techniques, using a synthesized data set. (a) shows the NPS calculated from data with no background trend 
(diamonds), and the same data but with linear and quadratic trends added (squares), (b) compares the data with linear and quadratic trends added before 
(squares) and after subtraction of a two-dimensional second-order polynomial fit to the data (triangles), (c) compares the data with linear and quadratic trends 
added before (squares) and after subtraction of a low-pass-filtered version of the data (crosses). 


n&x and wAy, corresponding to the reciprocal of the first 
zero crossings of the sine functions in the two dimensions, 
must be sufficiently large compared to the lowest frequency 
of interest in the NPS. In Fig. 3, S d w b is shown for n = m 
= 100 and for n = m = 200. Note that for m— 100, the first 
zero crossing of sinc(/wAyi;) is at a slightly higher fre¬ 
quency than the lowest frequency (u=l/(128A y)) in the 
spectrum. The effect is a slight reduction in the noise power 
at the lowest frequency, relative to that shown for m = 200. 
Because of this, in general, we have found the best results by 
making the size of the smoothing kernel somewhat larger 
than the size of the blocks used to calculate the individual 2D 
spectra. 

It should be noted that 2D quadratic fitting carries a sig¬ 
nificantly higher cost in processing time than does 2D 
smoothing. This is true even if the fits are performed on the 
relatively smaller blocks individually, rather than on the full 
data set. For example, on a Sun Sparcstation 10 in our labo¬ 
ratory, trend removal using 2D fits to 800 128X 128 blocks 
takes approximately 4 h of processing time, whereas only a 
few minutes are required for 2D smoothing of the same 
amount of data. 

C. Units and scaling 

The noise power spectrum is more precisely termed the 
power spectral density, that is, the noise power per unit 
frequency. 8 The quantity S(u,v)dudv is the contribution to 
the variance from noise components with spatial frequencies 
between u and u + du , and v and v + dv. S(u,v) thus has 


dimensions of inverse spatial frequency squared. The correct 
scaling of the NPS can be verified by observing that the 
volume (area) under the 2D (ID) NPS equals the total vari¬ 
ance. Thus, integration over frequency of spectra estimated 
using mean-zero sections from uniform illumination images 


Trend Subtraction 



Spatial Frequency (mm' 1 ) 

Fig. 3. Comparison of background trend removal techniques, applied to 
images obtained with one of the digital mammography systems, obtained 
using nominally uniform irradiation. The graph shows four ID power spec¬ 
tra, S( 0,u), obtained from thin cuts from four different 2D spectra. The 2D 
spectra were calculated following (a) no trend removal (none), (b) subtrac¬ 
tion of a fitted quadratic surface (fit), (c) subtraction of a low-pass-filtered 
version of the data using a 100X 100 pixel convolution filter, and (d) same 
as (c), but using a 200X200 pixel filter. 
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should numerically equal the mean-square deviations from 
their average of the pixel values in those sections. 

Noise power spectra are sometimes presented normalized 
by the square of the large-area (average) output signal. 20,21 
This normalization has the effect of compensating for differ¬ 
ence in gain between two systems which are being com¬ 
pared, and of simplifying the expression for the DQE. With 
this normalization, Eq. (6) becomes 


DEQ(w,u,<p) = 


C(w,u)MTF 2 (w,i;) 

5(« )t ;,<p)SNR? n 


Cju, i>) 
SNR? n 


NE Q(u,v,<p) 9 


(27) 


where S(u,v,(p) = S(u,v,<p)/(<pg) 2 is the normalized NPS, 
and NEQ(w,i?,<p) is the noise equivalent quanta. 


IV. RESULTS 

A. Two-dimensional NPS 

Neglecting digitization noise, the additive noise in the 
NPS obtained from a CCD-based digital mammography de¬ 
tector can be written as 

cr 2 a = a 2 + cr 2 r + S p (u y v ), (28) 

where a t is the CCD thermal (Johnson) noise, a r is the read 
noise, and S p (u,v) is the variance due to frequency- 
dependent periodic noise. These noise contributions are 
summed in quadrature because they are independent with 
respect to each other. Periodic noise can arise from capaci- 
tively coupled clock pulses, electromagnetic interference 
from power supplies, or 60 Hz pickup due to ground loops. 
The location in the image of the noise due to such periodic 
interference is not fixed from frame to frame, because its 
phase (in time) is not fixed with respect to detector readout. 
It is, therefore, not removed during the uniformly correction. 
Note that spikes from periodic additive noise are apparent on 
the v axis (the central, vertical axis in each of the plots in 
Fig. 4), but not on the u axis. As is shown below, these same 
structures are apparent in the Id spectra obtained by scan¬ 
ning a synthetic slit in the y direction. 

2D spectra obtained from both developmental systems ex¬ 
hibit strong periodic structures over all exposure ranges 
tested. These structures take the form of sharp peaks in the 
spectra, and are due to periodic additive noise. Figure 4 
shows gray-scale plots of the 2D spectra of one digital mam- 
mographic system at four exposure levels. Gridlines are 
spaced at 1 mm -1 . With increasing x-ray exposure, the peri¬ 
odic noise peaks are decreasingly apparent as their magni¬ 
tude becomes progressively less relative to that of the x-ray 
quantum noise contribution at low frequency. 

B. Effect of varying slit length on ID NPS 

We have studied the effect of varying the length of the 
synthetic scanning slit used for the estimation of S(w,0) or 
S( 0,u). Figure 5 shows ID spectra obtained from uniform 
illumination images using an exposure of 38 mR at the de¬ 
tector surface. Each of the spectra shown in the plot was 
generated by stepping a synthetic slit with x and y dimen¬ 



Fig. 4. (a) 2D NPS at 9 mR exposure to the detector surface. Gridlines are 
spaced by 1 mm -1 (b) 2D NPS at 18 mR exposure to the detector surface. 
Gridlines are spaced by 1 mm -1 , (c) 2D NPS at 38 mR exposure to the 
detector surface. Gridlines are spaced by 1 mm" 1 , (d) 2D NPS at 76 mR 
exposure to the detector surface. Gridlines are spaced by 1 mm" 1 . 


sions of 1 X L pixels, over the data in the x direction in steps 
of one pixel. The longest slit length shown (I = 512 pixels) 
corresponds to a slit approximately 20 mm long. At each 
step, the pixels in the slit were averaged. The resulting 128 
element series was then Fourier transformed and scaled as 
follows to give the measured ID NPS. Combining Eqs. (15) 
and (16) we have 


Scanned Slit, 38 mR 



Fig. 5. ID noise power spectra obtained from the same data set, but with 
various scanning slit lengths. Spectral components above approximately 0.6 
cycles per millimeter are insensitive to changes in slit length for slit lengths 
of 64 pixels or greater. However, a slit length of 4 pixels results in under¬ 
estimation of the noise power out to ~3 cycles per millimeter. 
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(29) 


where p'(n Ax) = 1/AT^ 2 ^i 0 1 p(nAxJby) is the average 

of the pixel values in the synthetic slit. Note that from Eq. 
(17), a section of the digital 2D NPS along the u axis, cal¬ 
culated from a NXM section of the image, is 

S(*AK,0) = lim WAx ^.^f Ay 


X 


N-\ 


M- 1 


2 


2 e - 2ni[{k * uHnAx)] 2 />(«Ax,mAy) 

n =0 m=0 


(30) 


For the special case when the slit length N L Ay is equal to the 
scan length WAx, Eq. (29) and (30) reduce to the same ex¬ 
pression: 


Ax Ay 

S(*Ah, 0)= lim — — 

N L Ay-*<x> 


X 


N-\ 


2 


2 e _2lrit( * A " )( " Ajc) W I p'(nAx) 

n — 0 


(31) 


Thus, for example, a cut along the u axis from a 128 
X 1282D NPS is identical to the ID NPS obtained by scan¬ 
ning a 1 X 128 pixel synthetic slit over 128 pixels in the x 
direction. In all cases, we verified that the ID NPS obtained 
from the scanned slit method with the slit length equal to the 
size of the data block used to generate the 2D NPS, was 
identically equal to the section from the 2D spectrum along 
the u or v axis. 

As shown in Fig. 5, the very low-frequency values of the 
1D NPS did not converge to a plateau value for slit lengths 
of up to 512 pixels, which is four times the scanned distance. 
This behavior is similar to that observed by Sandrik and 
Wagner for the Kodak Hi-Plus/XRP screen-film system. 18 In 
that study, measured NPS values at 0.39 cycles/mm contin¬ 
ued to increase with increasing slit length up to Z, = 5.9mm 
for scan a distance of 2.5 mm, while convergence to a pla¬ 
teau value was observed for other screen-film systems at 
similar frequencies. 


C. Slices from the 2D NPS 

Notably, in the context of the evaluation of storage phos¬ 
phor systems, ID spectral estimates have been obtained by 
extracting thick cuts from the 2D NPS. The cuts are made 
parallel to the primary axis of interest, but do not include the 
axis, so as to avoid low-frequency trending effects. We have 
evaluated this approach for analysis of the spectra of the 
digital mammography systems. Figure 6 is a plot comparing 
a thin cut (a single-frequency-bin wide) along the primary 


Cuts Along V-Axis, 76 mR 



Spatial Frequency (mm* 1 ) 


Fig. 6. Slices from 2D NPS of digital mammographic system. The solid line 
is a slice one-bin wide, along the v axis. The dotted line is the sum of two 
four-bin sections taken from either side of the v axis. Differences between 
the two occur primarily at low frequency and because of an-axis periodic 
noise components. 


( v ) axis of the 2D NPS with a thick cut generated by sum¬ 
ming the two four-bin-wide sections immediately on either 
side of, and adjacent to, the v axis. The thick slices ad¬ 
equately represent the on-axis slice except at low frequency, 
where the NPS decreases rapidly with increasing u. Of 
greater significance is the omission of the periodic noise 
structures near spatial frequencies of 3 and 6 cycles/mm. 
These represent noise power from fluctuations in the y direc¬ 
tion (left to right direction of the breast support) that are 
constant in the anterior-chest direction. Such noise may ap¬ 
pear as horizontal stripes in the mammographic image. 

As would be expected, thicker slices are poorer approxi¬ 
mations to the on-axis cuts, as the falloff in the 2D NPS in 
the dimension orthogonal to the slices causes the average 
spectral values to be reduced. This produces an effect similar 
to that resulting from utilization of a scanning slit of insuf¬ 
ficient length. 


D. NPS of the developmental mammographic 
systems 

Spectral estimates were obtained from images acquired 
using two developmental digital mammography systems lo¬ 
cated at the University of Virginia. Both systems contain 
butted arrays of CCD-based modules, but they differ in the 
type of x-ray converter and in the details of the optical cou¬ 
pling to the CCDs. The pixel size is ~37 pm for system 1 
and ~47 pm for system 2. The NPS of uniform illumination 
images from both digital mammographic systems was mea¬ 
sured over a range of detector entrance surface exposures 
from 9 to 78 mR. For each exposure, 2D noise power spectra 
estimates were made as described in Sec. Ill A. Figures 7(a) 
and 7(b) are semilog plots of cuts along the u and v axes of 
both systems at four exposure levels. The spectra have been 
normalized by the respective large-area output signal 
squared. Immediately apparent from Fig. 7 is the fact that the 
noise behavior of system 1 is different in the x and y dimen¬ 
sions, with the y dimension exhibiting significantly greater 
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System 1 



System 2 



Spatial Frequency (mm 1 ) 


Fig. 7. Noise power spectra of the two developmental digital mammography systems studied. Spectra are shown at four exposures for each system. In the 
spectra of system 1, the magnitude of 5(0,v) exceeds that of 5(w,0) at high frequency for all exposures, even though their magnitudes at low frequency are 
equal. 


noise at high frequency. In both systems, the y direction 
corresponds to the parallel transfer direction of the CCDs. 
Also, while the magnitude of the system 2 NPS decreases 
monotonically with increasing exposure, that of system 1 is 
higher at 76 mR than at 18 or 36 mR. This unexpected in¬ 
crease in noise is caused by the appearance of uncorrected 
fixed pattern noise in the images of system 1 at higher expo¬ 
sures. 

Also apparent are the very different shapes of the spectra 
obtained from the two systems, with the spectra of system 1 
leveling off at high frequency, while those of system 2 con¬ 
tinue decreasing out to the Nyquist frequency. At all expo¬ 
sures tested, the normalized noise power of system 2 is less 
than 10” 6 at Nyquist. As is shown below, attenuation of the 
high-frequency noise in system 2 is primarily a result of the 
smoothing effect of the spatial dewarping algorithm used. 


E. Comparison with pixel variance 

Each of the 2D spectra presented in this study was inte¬ 
grated over spatial frequency in order to compare the result¬ 
ing value with the variance calculated from the pixel-to-pixel 
fluctuations. The integrals were calculated over the range 
from —f N to +/w, where f N is the Nyquist frequency. In all 
cases, differences between the integrals and the directly cal¬ 
culated total variance were less than 0.5%. 

Figure 8 shows the total noise of the two digital mam¬ 
mography detectors over the range of exposures from 9 to 78 
mR. The noise was determined from the square roots of the 
integrated volumes under the 2D spectra whose axial slices 
are shown in Fig. 7. Figure 8 plots the log of the noise versus 
the log of the exposure. For both systems, straight lines have 
been fit to the experimentally determined data points, and the 
calculated slopes of the lines are shown in Fig. 8. In a system 
for which x-ray quantum noise is the dominant noise source, 
the slope of such a plot is approximately 0.5. A slope less 
than 0.5 indicates the presence of significant contributions to 
the noise power from sources other than x-ray noise. 
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F. Effects of image correction 

All detectors exhibit imaging characteristics that are im¬ 
perfect but constant from image to image. This type of im¬ 
perfection, including the presence of tumed-on or dead pix¬ 
els, local or global spatial distortions, and spatially 
fluctuating sensitivity, can be corrected in software. How¬ 
ever, correction algorithms can have an appreciable effect on 
the noise power spectra, and care must be taken in this regard 
when interpreting the spectra. Figure 9 shows two pairs of 
spectral estimates calculated from the same set of uniform 
x-ray exposures. In one pair, the raw data (prior to any cor¬ 
rections) were used. In the other case, the data were pro¬ 
cessed using the sequence of calibration steps typically ap¬ 
plied to clinical data: masking of bad pixels, subtraction of a 
dark frame, geometric dewarping, and sensitivity nonunifor¬ 
mity correction. For both raw and processed data, spectral 
estimates were calculated from difference images obtained 
by subtracting one image from another, pixel by pixel. The 
most apparent difference between the spectra of the raw and 
calibrated images is at high frequency, where the calibrated 


RMS Noise vs Exposure 



Log Exposure (mR) 


♦ System 1 
■ System 2 


Fig. 8. Log noise vs log exposure for two digital mammographic detectors. 
The noise was found from the square root of the integrated volume under the 
2D NPS. The plot shows the slope of straight lines fit to the measured data. 




















1290 


Williams, Mangiafico, and Simoni: Noise power spectra of images 


1290 


S(u,0) 72 mR 



Fig. 9. Comparison of the NPS obtained from corrected vs raw image data. 
The reduction in the high-frequency noise in the corrected image spectra 
occurs during dewarping. Of course, the detector MTF exhibits a similar 
reduction at high frequency. 

images exhibit much less noise power. The cause of the at¬ 
tenuation of high-frequency noise is the dewarping algo¬ 
rithm. During dewarping, the signal in the raw image is re¬ 
assigned to new locations in the dewarped image based on a 
matrix of displacement vectors stored in a look-up table. In 
general, a signal from any given pixel in the warped image is 
mapped to a location that does not coincide with the center 
of a pixel in the dewarped image. One method of handling 
this is to perform a proximity-weighted distribution of the 
signal among pixels whose boundaries fall within a half pixel 
of the corrected location. This has the effect of smearing out 
the signal content of each pixel in the raw image over several 
pixels in the dewarped image, a smoothing operation. 

Note that this method of dewarping is only one of many 
potential approaches, each of which can have a different ef¬ 
fect on the final image. An important function of noise power 
spectral analysis is to identify the impact of various image- 
correction algorithms on the image noise. For this reason, it 
is useful to obtain noise power spectral estimates both from 
corrected data and, when available, from raw data. 

G. Effects of aliasing 

Aliasing of noise components present above the Nyquist 
frequency is unavoidable because the only estimate of the 
NPS available from a digital system is that obtained follow¬ 
ing sampling by the digital matrix. In order to estimate the 
effect of aliasing in the measured NPS of the digital mam- 
mographic systems, cascaded linear systems analysis was 
used to calculate the presampling NPS. 6,29-33 Details of the 
theoretical development are presented elsewhere. 34 For both 
systems, the presampling NPS can be written as 

t- £ I mtf («^)I 2 + [(G4 > 0 ) + ^ 2 ] 

V^s 

X | sine (p x u ) sine (p y v) \ 2 + <r 2 , (32) 

where G is the system gain in CCD electrons per incident 
x-ray photon; d> 0 is the incident x-ray fluence per pixel; rj is 


Normalized NPS at 9 mR 



Fig. 10. Comparison of measured and theoretical digital NPS. To calculate 
the theoretical digital NPS, noise power in the theoretical presampling NPS 
above the Nyquist frequency was aliased into the frequency range between 
zero and Nyquist. The presampling NPS was calculated using cascaded 
linear systems analysis. Although the system gain, phosphor absorption ef¬ 
ficiency, and Swank factor are all energy dependent, their spectrally 
weighted values are approximately: G~3.5 ^-/incident x ray, 17 —0.50, and 
A- 0.90. 

the x-ray absorption efficiency; A s is the Swank factor; 
MTF(w,u) is the overall detector MTF; ctj is the rms CCD 
thermal noise in electrons; p x and p y are the pixel sizes in 
the x and y dimensions, respectively; and cr r is the rms read 
noise in electrons. G, V » and^ 5 are each functions of the 

x-ray photon energy, and bars above symbols denote 
weighted averages over the x-ray spectrum incident at the 
detector surface. Figure 10 shows an example of the aliased 
noise in images obtained at 9 mR exposure using one of the 
developmental systems. Plotted are the measured digital ID 
spectra, S^(w,0) and 5^(0,u) obtained from differences of 
raw images; the calculated presampling spectrum, 5^; and 
the calculated digital NPS, S digii ^ • In calculating the digital 
NPS, 5 pre values between one and two times the Nyquist 
frequency were added to those in the interval between 0 and 
the Nyquist frequency, according to the formula 

•^digital! w ) = 5 pre( w ) + 5 pre( 2a) N ~ 0>)~ a], 0< 0 

(33) 

The variance due to read noise is not aliased because it is 
added after the digital sampling process, and is merely an 
offset in the measured NPS. The calculated noise power 
above twice the Nyquist frequency is negligible. The close 
agreement in Fig. 10 between the measured NPS and the 
calculated NPS with aliasing shows that aliased noise can 
constitute a significant fraction of the high-frequency noise 
power in the mammographic systems. Note, however, that 
dewarping (c.f. Fig. 9) smoothes the high-frequency noise, 
thereby greatly reducing the effects of aliasing. 

V. DISCUSSION 

Similarly to other types of radiographic systems, digital 
detectors for mammography can usefully be characterized 
through Fourier analysis of their noise characteristics. How¬ 
ever, investigators preparing to estimate noise power spectra 
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of such system are confronted with a bewildering array of 
procedural considerations, lack of attention to which can lead 
to misleading results. Limitations formerly imposed by the 
tradeoff between spectral resolution and uncertainty in the 
spectral estimates are much less stringent than even a decade 
ago, since significant digital processing power is now avail¬ 
able to the majority of medical physicists. Thus, ensemble 
averaging over a large number of images obtained using uni¬ 
form x-ray exposure is feasible, thereby permitting fairly 
subtle features of the NPS to be resolved. The necessary 
resolution is determined by the proximity in Fourier space of 
relevant features and should be determined on a system-by- 
system basis. For the digital systems characterized in this 
paper, 128X 128 2D spectra, yielding 64 frequency bins in 
the single-sided ID spectra, were sufficient to resolve the 
features of interest. Occasionally, we found it useful to re¬ 
calculate particular spectra at higher resolutions. This may be 
useful, for example, for assurance that low-frequency struc¬ 
tures are, in fact, independent of background trend noise 
power residing in the lowest-frequency bin (the range be¬ 
tween zero and 1/ATAx, where N&x is the length of the data 
series). 

The noise in images from both the digital mammographic 
systems has markedly nonisoplanatic properties, particularly 
at low exposure where frequency-dependent system noise is 
most predominant. Such noise can be due to digital clock 
pulses that are radiatively or capacitively coupled to the ana¬ 
log data signals. Furthermore, the majority of this noise 
power resides away from the principle spatial frequency 
axes, and thus, does not appear on the 1D NPS, whether it is 
obtained from an extracted slice of the 2D NPS, or from a 
scanned slit. It is, therefore, important that the full 2D NPS 
be estimated, even if ID spectra are subsequently used for 
SNR calculations. 

The scanned slit method for obtaining ID spectra yielded 
results identical to those obtained from thin cuts through the 
2D spectra, except for the susceptibility of the former ap¬ 
proach to slit-length-dependent variability in the low- 
frequency components. The absence of convergence of the 
low-frequency values of S(w,0) reflects the peaked nature in 
the v direction of the 2D NPS, S(u,v), near v = 0. Since the 
value of S ss (w) obtained from scanning is scaled by L to get 
S(ufi) [c.f. Eq. (15)], only for values of u where the addition 
to the integral of Eq. (14) is proportional to L will the result 
be insensitive to changes in L. Thus, approximating the func¬ 
tion sinc 2 (Z,i>) as a triangle of base width 2/1 and unity 
height centered at u = 0, S(w,0) obtained using a slit of 
length Lj is approximately equal to that obtained with slit 
length Z, 2 only if S(w,u) is approximately constant in the v 
direction between v-2!L X and v = 2/L 2 - Note that except 
for w = 0, where S(0,u) is symmetric about the u axis, this 
requirement must also be met between u = — 21 L x and v 
= -2/Z, 2 . Dainty and Shaw have recommended that the slit 
length be at least equal to the inverse of the frequency reso¬ 
lution (that is, at least equal to the scan distance). 17 Our 
results indicate that for the digital mammography systems, a 
length significantly greater than the scan distance may be 
required for asymptotic convergence of the low-frequency 
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NPS estimates. In practice, the most direct method of deter¬ 
mining sufficiently large L is probably inspection of cross- 
sectional cuts through the 2D NPS, obtained at as high a 
frequency resolution as is practicable, bearing in mind the 
guidelines of the preceding sentences. 

For the noise data studied here, subtraction of a duplicate 
data set that has been convolved with a two-dimensional rect 
function kernel is an effective and efficient method for the 
removal of low-frequency trending effects. The effectiveness 
of the technique is indicated by the similarity between thick 
slices extracted from the 2D NPS near the central axes, and 
on-axis slices (c.f. Fig. 6). ID spectra obtained from thick 
slices parallel to the primary axis of interest may not be 
representative of the on-axis NPS, because the 2D NPS can 
contain on-axis features not present in adjacent regions of 
frequency space. For this reason, thin on-axis cuts are rec¬ 
ommended. 

Algorithms for correction of local detector imperfections, 
geometric distortions, or sensitivity nonuniformity can either 
amplify or attenuate components of the NPS. The best and 
most straightforward way to quantify the noise prior to post¬ 
acquisition processing is to obtain noise power spectral esti¬ 
mates from the difference of raw images. Here, raw means as 
soon as possible after analog-to-digital conversion. When 
evaluating commercial systems, it may not be possible to 
obtain data prior to image corrections without assistance 
from the manufacturer. 

It is our opinion that noise power spectral estimates 
should also be made with images processed using the soft¬ 
ware correction steps that are to be used in routine practice. 
Such data are always readily available, and provide the most 
realistic assessment of noise components present in clinical 
images. This analysis involves simply calculating spectral 
estimates from a number of individual processed images, and 
averaging spectra to reduce the uncertainty in the final esti¬ 
mate. In order to distinguish spectral components that arise 
from fixed pattern noise that is not fully corrected, these 
spectra should be compared to those obtained from the same 
data set, but subtracted from each other pairwise. The latter 
spectra contain no contributions from fixed pattern noise, 
since that noise is subtracted out (see Sec. IID 2). 

Calculation of an expression for the presampling NPS via 
cascaded linear systems analysis indicates that the shape of 
the 1D NPS is governed by three terms. The first term arises 
from x-ray quantum noise and the phosphor conversion sta¬ 
tistics, and is modulated by the system MTF. The second 
term arises from the Poisson noise associated with the CCD 
electrons generated by visible photons from the phosphor 
(secondary quantum noise) and generated thermally. This 
Poisson noise is modulated only by the pixel sampling aper¬ 
ture function. The third term is the frequency-independent 
(white) read noise. Comparison of the measured ID digital 
spectra with the theoretical presampling spectra shows good 
agreement for spatial frequencies up to ~6 cycles/mm. At 
higher frequencies, the digital NPS is increasingly larger 
than the theoretical presampling NPS with increasing fre¬ 
quency. This discrepancy is consistent with aliasing of the 
NPS components above the Nyquist frequency back into 
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those between zero frequency and the Nyquist frequency. 
The aliased noise consists primarily of secondary quantum 
noise, because the x-ray quantum noise is greatly attenuated 
by the system MTF above the Nyquist frequency. 

VI. CONCLUSIONS 

Fourier analysis has been used to characterize the noise 
properties of two developmental detectors for digital mam¬ 
mography. The noise power spectra of images from the two 
systems exhibit several common features. They contain fea¬ 
tures that are strongly dimension specific, and are often lo¬ 
calized to relatively small portions of frequency space. They 
also contain predominant peaks indicative of periodic addi¬ 
tive noise, both on and off the primary frequency axes. On 
the other hand, the noise properties of the two mammo- 
graphic detectors differ significantly in some aspects. Over a 
range of input exposures, the spectra of system 1 are charac¬ 
terized by nearly constant noise power for frequencies 
greater than about 8 cycles/mm, while those of system 2 
exhibit increasingly sharp falloff with increasing frequency. 
Cuts along the two primary axes of the 2D NPS of system 1 
exhibit very different shapes, while those of system 2 are 
nearly identical along the two axes. At all frequencies, the 
relative noise of system 2 decreases monotonically with in¬ 
creasing x-ray exposure, while that of system 1 reaches a 
minimum then begins to increase again. The overall noise of 
system 2 increases approximately as the square root of the 
number of x-ray quanta, indicative of x-ray quantum-limited 
behavior, while that of system 1 increases more slowly, in¬ 
dicating that a significant portion of the noise power arises 
from sources other than x-ray noise. 

The complex nature of the frequency distribution of the 
noise power in images from digital mammographic systems 
makes analysis of the complete 2D spectra, rather than only 
the 1D spectra, important. Straightforward removal of back¬ 
ground trending artifacts is possible via a simple filtering 
process. Cuts obtained along the primary axes of the 2D 
spectra can provide valid estimates of the ID NPS without 
the use of synthetic scanned slit techniques. Spectral esti¬ 
mates from difference images and averaged images can be 
used to isolate stochastic, fixed pattern, and varying nonsto¬ 
chastic components for identification of noise sources. 

Aliasing can contribute significantly to noise power near 
the Nyquist frequency. Cascaded linear systems analysis is a 
useful tool for evaluation of the magnitude and sources of 
aliased noise. 
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